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Achioving  balance  Lo  a data  aet  is  a desinable  feature  since  eetimacote  and  teats 
concerning  the  variance  coinponente  of  a gjven  model  have  optimal  properties,  par* 
ticuiarly  under  the  uormallty  assumptloD.  Complete  data  balance,  however,  ie  not 
always  attainable.  Therefore,  the  dfect  of  imbalance  of  a design  on  a particular  quan- 
tity of  interest  is  always  a concern  to  an  experimenter.  In  this  dissertation,  we  study 
the  overall  effect  of  data  imbalance,  as  well  as  the  effect  of  variance  components,  usng 
a methodfor  generating  designs  having  specifled  degrees  of  imbalance.  The  probabil- 
ity of  getting  a negative  ANOVA  estimate  of  a variance  component,  the  power  of  an 
F-test  concerning  a variance  component,  with  or  without  homogeneity  assumption 
of  error  variances,  and  the  coverage  probability  of  a confidence  interval  on  a variance 
component  are  considered  as  quancitiee  of  interest.  Generalized  linear  models  tech- 
niques are  used  to  model  a quantity  of  interest  in  terms  of  factors  considered  to  have 
influence  on  such  s quantity.  This  modeling  scheme  allows  us  to  acquire  a bettor 
iusight  into  the  effect  of  imbalance  on  the  quantity  of  interest  under  consideration. 


CHAPTER  1 
INTRODUCTION 


In  3D  experlmeDtal  situAtioa,  the  levels  of  a factor  not  be  subject  to  the  control 
of  the  eiporimonlor,  or  they  mto'  be  selected  at  candoin.  The  effect  of  such  a factor 
is  considered  random,  and  ite  variance  is  termed  variance  component.  Furthermore, 
data  sets  obtained  in  an  experiment  may  be  unbalanced.  By  this  it  is  meant  that  the 
numbera  of  observations  in  the  various  subclasses  of  the  data  ace  not  all  equal.  The 
pattern  that  these  numbeis  take  in  a given  experimental  Bituatlon  Is  referred  to  as  a 
design.  In  particular,  if  these  numbers  are  equal,  then  the  data  set  is  called  balanced, 
and  the  associated  design  is  said  to  be  completely  balanced.  The  design  used  in 
connection  with  a ^ven  model  plays  an  important  role  in  assessing  the  signiffcance 
of  the  various  effects  in  the  model.  This  is  particularly  true  in  the  analysis  of  the 
model's  variance  components. 

Several  designs  are  often  possible  for  a given  cost  of  experimentation.  Determina* 
tion  of  optimal  designs  and  estimation  of  parameter  of  a fitted  model  are  therefore 
needed.  The  choice  of  a deaign  and  the  corresponding  method  of  analysis  normally 
depend  on  which  functions  of  variance  components  are  estimated.  While  much  has 
been  done  in  the  area  of  estimation  of  variance  components,  design  aspects  associ- 
ated with  the  analysis  of  variance  components  are  less  developed.  This  discrepancy 
arises  in  part  from  the  lack  of  a general  methodology  to  compare  designs.  In  sev- 
eral studies,  designs  were  compared  using  various  criteria-  Meat  of  these  criteria, 


r.  employed  a single-valued  Function  such  as  the  trace  or  determinant  oI  the 


Achic%ing  balance  In  a data  net  Is  a desirable  feaciire  due  to  Che  opcimal  prop- 
erties associated  with  estlmatois  and  tests  concerning  the  variance  components  of  a 
given  model  under  ihe  normality  aseumptlon.  Complete  data  balance,  however,  is 
not  always  attainable.  Therefore,  there  is  a need  to  have  at  one's  dispceal  adeqnaCe 
procedures  tar  the  analysis  of  nnbalanced  data.  The  analysts  of  unbalanced  models 
is  quire  eomplicaled.  Unlike  balanced  models,  there  is  no  unique  analysis  of  variance 
(ANOVA)  sini»  the  partitioning  of  the  total  sum  of  squares  can  be  done  In  a variety 
of  wuys.  Moreover,  the  distributional  properties  of  chi-squaiedness  and  independence 
of  the  sums  of  squares  are  no  longer  valid.  As  a result,  rrptimal  properties  of  the 
ANOVA  estimaconi,  namely,  nrinunuju  variance  quadratic  unbiasedness  and  uniform 
minimum  variance  unbiasedness  under  normality  are  lost  (see  Seaile,  Casella,  and  Mc- 
Culloch, 1902,  pp.  221-225).  Therefore,  whenever  we  have  an  unbaianrsd  data  set, 
questions  arise  os  to  how  unbairuiced  the  data  set  is,  what  effects  does  imbalance  (rio- 
balancedness]  have  on  data  analysis.  To  answer  the  first  question,  we  need  a certain 
measure  for  determining  the  degree  of  departure  from  balance.  ^^^lCh  this  measure,  we 
can  quantity  the  degree  of  imbalance  rather  chan  using  subjective  descriptions  such 
as  “irxlremely  imiNtlanced,"  “moderately  unbalanced,"  or  “nearly  balanced.”  Khuri 
(1987.  p.  395)  Slated  that  such  a measuni  can  serve  “as  an  indicator  of  the  suitability 
of  approximate  metliods  that  are  adapted  from  balanced-data-based  procedures  and 
used  to  anatitfie  an  unbalanced  model." 

Unbalanced  designs  present  dlfiiculties  to  comparative  studies  of  methods  for  esti- 
mating variance  compooenu.  Even  the  simple  Idea  of  comparing  sampling  variances 
is  intractable  because  the  variances  themselves  are  not  only  fluictinns  of  the  variance 
components  but  also  functions  of  the  design,  D (Searle,  Caselia,  and  McCulloch,  1992, 
p.  222).  In  this  context,  the  main  dilficulty  in  deeigning  studies  Is  to  “adequately  cover 
all  possible  designs"  so  as  to  provide  a basis  for  drawing  general  conclurnons.  One 
may  choose  a collection  of  patterns  of  cell  frequendes  “which  Intuitively  seem  to  range 


from  sli|ht]y  lo  badly  unbalaoced"  (Swallow  and  Searle,  197S^  p.  270).  However,  as 
will  be  been  later,  having  a method  for  measuring  the  degree  of  Imbalance,  coupled 
with  the  ability  to  generate  designs  with  a specified  degree  of  imbalance,  will  enable 
us  to  adequately  study  the  effect  of  imbalance  on  particular  quantities  of  interest.  For 

cance  arc  commonly  of  internst.  In  addition,  we  can  also  perform  better  comparisons 

Mathew,  and  Siuha,  1998,  Ch.  3,  pp,  58*88). 

In  this  dissertation,  we  have  concentrared  on  comparing  designs  for  estimating 
variance  components  using  a graphical  approach,  the  so-called  quantile  disperaion 
graphs  (QDGs)  which  was  first  introduced  in  Khnri  (1997).  We  have  also  investigated 
the  effect  of  imbalance  on  various  quantities  of  Interest  which  are  modeled  usmg 
generalized  linear  models  fechniques.  For  this,  the  method  of  generating  designs  with 
a specific  degree  of  imbalance,  as  introduced  in  Khuri  (1996).  is  used. 

This  dissertation  Is  divided  into  five  chaplets.  Chapter  1 fpves  a brief  iolroduc- 
tion  to  the  area  of  design  comparisons  and  the  study  of  the  effect  of  imbalance.  In 
Chapter  2,  we  review  the  pertinent  statistical  literature.  A related  survey  is  found 
in  Kliuri  and  Sahai  (1985.  pp.  298-294).  Chapter  3 is  devoted  to  comparing  desgns 
for  the  one-w^  random  model,  which  is  in  Section  3.2,  for  the  two-way  cress  classi- 
fication random  model,  which  is  in  Section  3.3,  and  for  the  three-fold  nested  random 
model,  which  is  in  Section  3.4,  using  the  QDGs  and  the  so-cailed  empirical  quantile 
dispersion  graphs  (EQDGs).  in  each  ease,  two  estimation  methods,  namely,  ANOVA 
and  ML,  are  considered  to  compare  these  estimation  methods  for  a given  deigo.  in 
Chapter  4,  an  investigation  is  made  concerning  the  effect  of  imbalance  on  certain 
quantities  of  interest,  including  the  probability  of  a negative  ANOVA  estimate  of  a 
variance  component  in  Section  4.1,  the  power  of  an  f-lest  for  a variance  component 


under  homogeneous  error  varionoes  in  Section  4.2,  the  power  of  an  F-ieot  for  n vari- 
ance component  under  heterogeneous  error  variances  in  Section  4.3,  and  the  coverage 
probahility  of  a eonfidenee  interval  for  a variance  component  b Section  4.4.  Finaliy. 
Chapter  5 presents  a discussion  of  future  research  issues.  The  program  for  generating 
designs  with  a specified  degree  of  imbaiance,  which  is  developed  for  this  research,  is 
provided  in  the  Appendix. 


CHAPTER  2 
LITERATURE  REVIEW 


2T  Dgsigp  Compariflons 


(I9S4),  whose  results  were  Isier  published  by  Anderson  and  Crump  (1967).  Anderson 
end  Crump  derived  on  optimel  nllocatioo  procedure  to  estimate  (iheamong-rlasses 
variance  component)  or  the  ratio  7}=~ 


Their  procedure  minimizs  the  variances  of  the  ANOVA  estimatois,  d\  and  ij,  for  a 
fixed  number  of  groups  (A)  and  a fixed  total  sample  »ee  (iV),  They  also  presented 
procedures  for  determining  the  optimal  choice  of  k to  estimate  and  p.  They  defined 
the  efficiency  of  estimating  or  p as  the  ratio  of  the  variance  of  the  estimate  with  an 
optimal  k (i.e.,  an  optimal  dmign)  to  that  of  the  estimate  at  a given  k (i.a,  a given 
de^).  They  showed  that  a suitable  choice  of  the  sampling  plan  (in  terms  of  A)  can 
provide  large  gains  in  the  cfliclencies  rrf  the  estimates,  and  p. 

Thompson  and  Anderson  (1975)  compared  several  balanced  deigns,  for  the  one< 
way  random  model,  using  ANOVA  estimators,  truncated  ANOVA  estimators,  max- 
imum lihelihood  (ML)  irstimaCora,  and  modilled  mrrximum  lihelihood  (MML)  esti- 
rnatOTS  of  the  variance  compnnenta.  The  last  estimators  were  intrrxluced  by  Klotz, 
Milton,  and  Znck.s  (1969).  Thomson  and  Anderson  also  considered  an  unbalanced 
design  which  allocates  one  observation  to  each  of  A|  classes  and  rwo  observations  to 
each  of  As(=A-Ai)  classes.  Such  a design  was  referred  to  as  a (A| , Az}-dcaign.  This 
design  was  found  to  be  optimal  for  estimating  when  p > L by  Andersoa  and 


Crump  (1967).  Thompson  and  Anderson  found  that,  regardlesn  of  how  large  ij  was, 
ML  provided  an  estimator  of  with  a smaller  asymptotic  variance  than  the  proce- 
dures of  obtaining  truncated  weighted  ANOVA  and  truncated  unweighted  ANOVA 
estimators.  However,  the  dilference  between  the  unw^ghted  ANOVA  and  the  ML 
results  was  small  when  rj  was  large. 

.Monte-Cario-based  comparisons  of  estimators  of  variance  components  with  un- 
balanced designs  were  mode  by  Swallow  and  Monahan  (1984),  and  Cbaloner  (1987). 
Both  papers  considered  one-way  random  model  with  an  arbitrarily  chosen  unbalanced 
(iesigiis.  Chaloner  showed  that  Bayesian  estimators  were  viable  competitoie  to  many 
clafiical  alternative  estimatots. 

Comparisons  of  designs  for  the  two-way  cross  clasDBeatioa  random  model  are  con- 
sidered  in  G^lor  (1060).  Substantial  contributions  In  this  area,  however,  were  made 

R.  L.  Anderson  and  his  colleagues.  Muse  and  Anderson  (1978)  compared  several 
particular  designs  using  ML  in  a model  without  interaction.  The  asymptotic  behav- 
ior of  ML  estimators  of  variance  components  was  examined  through  its  large  sample 
variance-covariance  results.  They  concluded  that  prior  iofonnalion  was  necessary  for 
selecting  a good  design  for  the  estimation  of  variance  components.  Using  the  trace 
of  an  asymptotic  variance-covariance  matrix  as  a criterion,  they  also  found  that  the 
balanced  design  was  best  only  when  of  is  the  dominant  variance  component,  that 
is,  when  of  is  larger  than  both  of  and  of,  the  variance  components  for  rows  and 
columns,  rest>ectively.  Conclusions  favoring  certain  unbaianced  designs  over  balanced 
designs  can  also  be  found  in  Mostafa  (1967).  Moslafa  considered  Che  two-way  croased 

Muse,  Anderson,  and  Thilehamol  (1982)  extended  the  result  of  Muse  and  An- 
derson (1978)  to  include  additional  comparisons  of  designs  for  the  two-way  cross 
ciassiheation  random  model  with  interaction  using  asymptotic  maximum  likelihood 


s.  Simil&r  findings  and  n«ulta, 


in  Muse  and  Anderson  (197S), 


Leone  end  Nelson  (1966)  employed  the  balanced  three-rold  nested  random  model. 
They  conducted  empirical  studies  of  the  sampling  distribution  ofANOVA  cstimalois 
using  eight  dilTerent  cumbinatious  of  variance  components  sets.  Histograms  of  6^ 
under  the  normality  assumption  were  provided.  They  also  showed  empirical  and  an* 
alytical  results  concerning  percentages  of  negative  estimates  of  variance  components. 

Leone,  Nelson,  Johnson,  and  Bisenstat  (1966)  considered  the  unbalanced  three* 
fold  nested  random  model.  They  employed  the  ‘^laggeted"  and  "inverted"  four-stage 
nested  designs  and  presented  empirical  results  concerning  the  distributions  of  the 
vnriance  components  for  both  detngns  with  eight  sets  of  variance  components.  The 
normal,  rectangular,  and  exponential  underlying  distributions  were  considered.  No 
particular  difference  between  staggered  and  inverted  deagns  was  detected  on  the  basis 
of  ths  emplricaJ  shape  of  the  distributions  of  the  .ANOVA  estimate  of  o^.  With  the 
normal  and  lectangular  underiying  distributions,  variance  components  estimates  from 
the  first  three  stages  behaved  quite  rimilarly  in  their  empirical  distribution.  However 
the  exponential  distribution  {iroduced  imprecise  variance  estimats  because  of  its 
tong  tail.  Since  Balubridge  (196S,  1965)  introduced  the  staggered  design,  it  has  been 
received  continuous  attention  by  several  researchers.  (See  Smith  and  Beverly.  1981; 
Nelson,  1993,  1999a,  1995b;  Khattrcc  and  Naik.  1995;  Khattree,  Naik,  and  Mason, 

Goldsmith  and  Gnylor  (1970)  considered  the  two-fold  nested  random  model  with 
variunce  comiionenrs,  and  of.  They  stutlied  the  effect  of  design  on  the 

precision  of  the  ANOV.A  estimates  of  these  variance  cumpooenta  using  computer 
simulatiOQ.  A total  of  61  unbalanced  designs  were  chosen,  and  for  each  demgn,  49 


were  considered.  Three  criteria  were  proposed  to  compare  designs  and  to  determine 


optimum  demgnsovEr  the  variance  component  configurations  and  sample  sizes.  These 
criteria  were  based  on  the  trace,  the  determinant,  and  the  weighted  trace  of  the 
vaiiance<covariance  matrix  of  the  vector  of  variance  component  cetimatora.  The  trace 
criterion  was  recommended  since  it  wss  sensitive  to  changes  in  the  sample  size  and  the 
variance  component  configuration,  as  well  as  being  the  easiest  criterion  to  compute. 
When  both  j}i(=  tnd  were  small,  a balanced  design  was  found  to  be 

optimal  in  tbe  sense  of  pioviding  tbe  smallest  value  of  the  trace.  When  rji  was  large, 
designs  with  the  highest  second*stage  degrees  of  &esdora  were  favorable;  when  tjg  was 
large,  designs  with  liigh  first-stage  sampling  were  seiKted. 

Cummings  and  Gaylor  (L974)  considered  variance  components  testing  for  an 
unbalanced  two-fold  nested  random  model.  An  exact  test  for  the  hypothec  of 
can  be  formed  regardless  of  unbalancedness  of  the  design,  but  not  for 
the  hypothesis  of  Ho:a^=0.  It  is  because,  with  an  unbalanced  design,  the  properties 
of  independence  and  chUsquaredness  of  mean  squorm  no  longer  hold  under  the  usual 
normality  assumptions.  Ignoring  this  failure  of  properties,  they  employed  Sattertb- 

teat  is  unavtulabic  In  balanced  AITOVA  table.  They  considered  "test  size  disturbance" 
as  the  discrepancy  of  the  estimated  size  of  this  approximate  test  from  the  assumed 
lest  size,  rv.  Dependence  and  non.chisquorednesa  of  a mean  square  arc  considered 
as  contributing  foctors  for  the  test  size  disturbance.  They  noticed  that,  in  the  case 
of  = n,y,  j jS  /,  j = 1,2,"*  ,k  [Khuri  (1997)  used  the  termiuoiogy  '‘partially 
balanced  nested  dc^gn"  for  this  case],  independence  of  the  mean  squares  Is  obtained, 
but  not  chi-squaredness.  In  this  ntnatlon,  the  test  size  disturbance  bereased  as  171 
increased.  On  the  other  hand,  when  designs  have  chi-squared,  but  dependent  mean 
squares,  the  test  disturbance  decreased  as  fh  increased  with  rebtlveiy  the  same 


mfignjtude.  When  both  properliea  of  independence  and  chi-squaredneaa 


laied,  the  test  siec  disturbance  was  apparently  small,  a bebavioc  which  Cummings 
and  Gaylor  (1974,  p.  769]  attributed  to  a so-called  "counterbalandog  effect." 

Tan  and  Cheng  (1984)  considered  the  problem  of  testing  a variance  component 
in  an  unbalanced  two-fold  nested  random  model.  They  employed  two  types  of  un- 
balanced designs  in  order  to  examine  and  compare  four  test  statistics  for 
with  respect  to  power  of  the  test.  They  obtained  approximate  power  functions  using 
Laguerre  polynomial  expan^ns. 

Mukerjee  and  Huda  (1988)  considered  raultifactor  random  effects  model  to  esti- 
mate variance  components  by  the  method  of  unweighted  squares  of  means.  They 
showed  that,  when  the  numbers  of  levels  of  the  tnelois  are  considered  fixed,  a bal- 
anced design,  if  it  exists,  is  optimal  in  the  sense  of  minimizing  variances  of  the  variance 

Khuri  (1997}  obtained  the  exact  distribution  of  an  ANOVA  estimator  of  a vari- 
ance component  by  deriving  its  quantiles  in  order  to  assess  the  quality  of  ANOVA 
estimation.  He  demonstrated  this  iraug  the  balanced  two-way  cross  dasmfication  ran- 
dom model  with  Interaction.  Since  the  expected  mean  squares  in  the  corresponding 
.ANOVA  table  are  functions  of  the  variance  components,  the  exact  distribution  of  a 
scaled  value  of  an  ANOVA  eetimator,  namely,  can  be  expresed  as  a linear 

combination  of  independent  chi-squared  variates,  and  is  a function  of  the  following 
ratios  of  variance  components:  8)=^,  02=^.  Using  the  algorithm  given  by  Davies 
(1980),  the  quantiles  of  H'.  namely  «(p),  that  satisfy  P[W  are  calculated 

for  each  p (p  e [0, 1|).  The  distribution  of  W can  be  depicted  umng  a quantile  plot 
of  tv  fbt  the  given  design  and  tbc  actual  values  of  0\  and  Qi-  Quantile  plots  provide 
information,  concerning  the  efficiency  of  by  examining  the  spread  in  the  values  of 
H',  the  extent  to  which  an  ANOV.A  estimator  can  be  negative,  and  the  stability  in 
tbc  values  of  IV.  Khuri  also  proposed  the  so-called  qusnlde  disptrsion  graphs  (QDGs) 
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not  only  as  a means  for  assessing  the  overall  quality  of  ANOVA  estimation  for  a par- 
ticular design,  but  also  as  an  effective  tool  to  compare  cbe  capability  of  competing 
designs  for  estimating  variance  components.  Maximum  and  minimum  quantiles  of 
U'  were  obtained  over  a set  of  scvernl  value  of  d~(ffi,d3j'  selected  in^e  a region  In 
tbe  parameter  space  of  6.  Tlie  QDGs  were  constructed  by  plotting  these  maximum 
and  minimum  quantOe  values  against  p.  Extension  of  this  methodology,  however,  is 
needed  to  other  estlmatom  of  a variance  component,  e.g.,  ML  and  REML  estimators. 

2,2  M«aSU""«  ris'e  lmh.l.ne« 


Low  (1976]  implicitly  mentioned  that  the  quantity  - it)*  could  measure 

devistioo  from  complete  balance  for  tbe  one-way  model  ^ce  it  would  be  small  for 


nearly  equal  rq’s-  Here,  rq  denotes  the  number  of  observations  from  the  i**  treat- 
ment. and  A = ^ fq.  Hess  (1979),  however,  was  the  Srst  to  propose  measuring 
imbalance,  thereby  extending  Low’s  idea,  in  his  study  of  the  sensitivity  of  minimum 
norm  quadratic  unbiased  estimatom  with  respect  to  tire  a priori  weights  for  <7^  and 
tr*.  the  model's  variance  components. 

Ahrens  and  Piocua  (1981)  made  the  first  formal  attempt  to  quanli^  the  degree 


of  imbalance,  developing  two  m 
number  of  groups): 


es  of  imbalance  for  the  one-way  model  [k  is 


where  D denotes  a design,  that  Is,  a set  of  rq-values  for  the  one-way  model.  The  first 
measure  is  such  that  0 < 'r(Ll)  < 1,  whereas  for  the  second  measure,  | < f{D)  < 1. 
Large  values  of  these  measures  indicate  *hiear  balance",  while  vnines  close  to  the 
respective  lower  bounds  indicate  severe  inibelance.  Ahrens  and  Pincus  utilized  these 


mefi^uRS  to  888068  the  efficiency  of  an  unbalanced  design  as  compared  to  a balanced 
design  nith  the  same  number  of  observations. 

Khurt  (1987)  proposed  a measure  of  data  Imbalance  as  a function  of  Pearson's 
chi'Squared  statiatic,  for  meaanring  departure  of  the  subclass  hequeocies  of  the 
dots  from  "complete  balance".  Ill  addition  to  complete  balance,  where  frequencies 
are  equal  within  all  Che  subclasses,  he  also  considered  certain  types  of  ''balance", 
such  as  "proportional  subclass  frequencies"  in  crosH-claaeilication  models,  "partial 
balance"  and  "last-stage  uniformity"  in  nested  models.  His  procedure  tnenKurcs  the 
degree  of  imlutlance  in  the  following  way:  based  on  an  appropriate  loglinear  model 
fot  the  subclass  (feqnencies,  ML  esiimatea  of  the  expected  frequencies  ace  calculated 
for  each  type  of  "balance" . Using  these  estimates.  Pearson's  approximate  chi-squared 
statistic,  which  is  used  in  the  goodness-of-fit  test  for  the  loglinear  model,  is  obtained. 
A sunple  function  of  that  statistic  leads  to  the  measure  of  data  imbalance,  which 
quantifies  departure  from  the  type  of  “balance"  considered.  For  the  one-way  model, 
the  measure,  t(D),  proposed  by  Khuri,  coincides  with  Abtens  and  Pincus'  (1981) 
nicaaure,  c(Z)).  The  generalixation  of  the  procedure  proposed  by  Kbuci  for  measuring 
data  imbalance  to  higher-order  models  is  straightforward. 

Lera  Marques  (1994)  extended  the  measures  of  imbalance  given  by  Ahrens  and 
Pincus  (1981)  to  higbec-order  models.  She  used  an  arithmetic  average  - wlthoui  pro- 
viding ony  rationale  for  this  cxtensloD  - of  Ahrens  and  Pincus'  measures  obtained  from 
each  stage  (for  nested  model),  or  from  each  row  and  column  (for  crossMlassification 
model).  She  lUuatrated  her  measures  uung  four  modals,  namely,  two-fold,  three- 
fold nested  models,  two-way  and  three-way  cross  clas»fication  models.  Since  an  un- 
weigiited  average  of  the  measures  of  imbalance  does  not  fully  convey  the  information 
of  data  structure.  Lera’s  procedure  should  be  modified  by  using  weighted  averages. 


23  Generation  of  Dacpn  md  mi  InvwtitHioii  of  tht  Effect  of  D«M  Imbtitnce 


To  study  the  effect  of  imbalance  on  data  analysis,  several  quantities  of  interest 
hate  been  considered.  The  following  is  a selected  Ust  of  these  quantities:  thecoverage 
probability  of  a confidence  interval  lor  a variance  component  (Thomas  and  Hultquist. 
1P7S);  the  ratio  of  mean  squared  errors  (MSE)  of  minimum  norm  quadratic  (MINQ) 
estimators  Ibr  completely  balanced  and  unbalanced  desigrra  (Abrens  and  Sanchez, 
19S3);  bias,  sampling  variance  and  hfSE  of  ANOVA  and  ML  estimators  of  variance 
componenta  (Caro,  Grossman  and  Fernando,  1905);  large  sample  relative  efficiencies 
as  the  ratios  of  variances  of  ANOVA  estimators  to  those  of  ML  estimators  of  rr^  and  of 
for  the  one-way  model  (Donner  and  Koval,  1907);  empirical  powers  of  the  F- test  and 
the  likelihood  tatio[LR)-test  of  of  in  the  one-way  random  model  (Donner  and  Kovrd, 
1989);  the  relative  efficiency  as  the  ratio  of  rha  variance  of  an  ANOVA  estimator  to 
that  of  minimum  norm  quadratic  unbiased  (MINQU)  estimator  of  rrf  (Ahrens  and 
Santrhez,  1992);  the  scaled  varirmce  of  an  ANOVA  estimator  of  of  (Khuri,  1998)  and 
the  tatirrs  of  empirical  powers  of  tests  for  of  (El-Basaiouni  and  Seely,  1996).  Except 
for  Khuri  (1996),  these  studies  were  based  on  certain  unbalanced  desigrrs  chosen  in  a 
subjective  manner  turd,  therefore,  lai)  to  provide  an  overall  coverage  of  other  possible 
unbalanced  situations.  Except  for  Caro,  Grossman  and  Fernando  (1985),  which  used 
a two-fold  nested  rnrtdom  model,  all  the  aforemeutioned  studies  employed  the  one-way 
randoiit  model. 

To  study  the  effect  of  imbalance  on  a quantity  of  interest,  a method  for  geoeratiog 
designs  having  s specific  degree  of  imbalance  is  needed.  Three  such  methods  were 
introduced.  Ahrens  aud  Sanchez  (1992)  conrndered  the  one-way  ratrdom  model  and 
irsed  one  of  two  measures  of  imbalance  proposed  by  Ahrens  and  Pincus  (1981)  to 
compare  de»gns.  Ahreits  and  Sanchez  uxambed  the  relative  eScieucy  of  ANOVA  es- 
titirators  with  respect  lo  MINQU  estimators  in  terms  of  the  ratio  of  variances  of  these 


cstiiiialore.  They  gmersled  designs  Iqr  simply  msuipulatliig  several  combinations  of 
sample  sizes  ranging  from  the  least  to  the  most  unbalemced  situations  for  a given  total 
number  of  oltservationn  {N)  and  a given  number  of  groups  (fe)  with  a fixed  first  group 
sample  size  (ni).  They  were  thus  able  to  compare  designs  with  diferent  degrees  of 
imbalance  after  fixing  N,k,  and  n,.  However,  they  were  unable  to  compare  "equally 
unbalanced"  designa,  i.e.,  desigus  with  the  saiue  depee  of  Imbalance  for  fixed  IV  and 

Donnet  and  Koval  (1587)  developed  a method  for  generating  designs,  UKog  the 
one-way  random  model,  based  on  an  underlying  probability  distribution,  such  as  the 
truncated  negative  binomial  distribution,  truncated  below  one.  They  Investigated  the 
large  sample  relative  eliiricncy  [RE)  of  AXOVA  estimators  ( jj,  e?)  with  respect  to  the 
,ML  estimators  based  on  their  variances  using  Monte-Carlo  simulatioo.  Note 

that  the  equation  for  RE{d^\d^)  iu  Donnet  and  Koval  (1987,  formula  3.1,  p.  184) 
waa  incorrectly  cited  from  that  of  Crump  (1551,  equation  on  p.  10).  Examining 
the  beliavior  of  ML  estimators  was  particularly  interesting  since  there  is  no  explicit 
form  of  tliese  csCimalors.  On  the  basis  of  a Mooie-Catlo  study,  they  reached  the 
conclusion  that  the  effect  of  imltalHiice  on  A£l(d^|5j)  depends  on  the  idioice  of  the 
intraclass  correlatiOD  coeificienc  (p=prfcs)-  As  a measure  of  imbalance,  they  used  one 
of  Ahrens  and  Pincua'  measures,  namely  d=i/(D).  One  of  the  main  deficiencies  of 
their  method  is  that  the  geiierauxl  design  docs  not  correspond  exactly  to  a specified 
degree  of  imbalance^  rather,  it  has  on  the  average  chat  degree. 

Khuri  (1956)  proposed  a new  method  for  constructing  designs  for  the  one  way 
model  having  a specified  degree  of  unbalance.  More  specifically,  pven  a measure  of 
imbalance,  e.g.,  d=v(D),  a design  is  to  be  generated  having  a specified  value  of  p 
and  a specified  value  of  rq)  for  a given  k.  Thus  the  design  D={ni, 

should  satisfy  n,=N  and  r^=^.  By  letting  s,=^,  0 < xt  < 1,  the  set 


of  &I1  TeA^ble  solutions  consists  of  the  coordinates  of  the  points  in  the  interior  of 
« regular  (*-  1 )-dimenaonal  simples  that  belong  to  the  inlerecction  of  a (t-  1)- 
dimensiona!  hyperplane  and  a (h—  l]<dimcnslonal  hypersphere.  This  intersection  is 
a (A  - 2)-dimertelonaI  hypeisphere  with  a radius  centered  at  the  centroid 
of  the  simples.  Since  a set  of  fessible  solutions  Is  s subset  of  this  hypcmphere,  the 
values  of  X,  having  a specified  value  of  d were  obtained  uaog  spherical  coordinates 
associated  with  the  hypersphere.  Once  values  of  X|,Xe,  • ,xs  are  determined,  if  NXi 
is  a positive  integer  for  i = 1,2, •*'  ,A,  then  = Nxt,  i = 1,2, • • • ,A.  produce  a 
design  having  the  specified  values  of  S and  d.  Otherwise,  the  equations  n,^N 
and  do  not  have  an  exact  integer  solution.  Khuri  showed  how  to  obtain 
an  approslinaCc  solution  to  these  equations  in  an  optimal  manner.  The  proposed 


design  would  be  {rt|,<  • • ,nj,}. 

As  a quantity  of  interest.  Khuri  (1996)  used  a scaled  variance  of  the  ANOVA 
estimator  of  oj  for  the  one-way  niodel,  namely,  the  quantity  g=jt  V'or(cJ)  in  order 


to  determine  the  effect  of  data  imbalance  on  the  precision  of  estimating  Values  of 
Q were  obtained  using  designs  generated  with  specified  combinations  of  N.  d,  and  p 
values.  An  empirical  relationship  was  established  between  ;,  on  one  hand,  and  JV.  d, 
and  p,  on  the  other  hand.  This  was  done  by  fitting  a response  surface  model  in  which 
g is  the  response  variable  and  Af,  d,  and  p are  the  control  vaxiables.  By  exploring 
the  response  surface  corresponding  to  this  empirical  relationship  over  the  factor  space 
associated  with  the  specified  combinalioiis  of  Al,  d,  and  p,  it  is  possible  to  assess  the 
effects  of  these  parameters  on  the  precision  of  estimating  o^. 


power  of  tests.  Donner  and  Koval  (1969)  investigated  the  effect  of  imbalance  on 
the  power  of  the  T-test  for  the  variance  components  in  the  one-way  random  model. 
Using  the  technique  piopceed  by  Donner  and  Koval  (1987)  to  generate  designs  for  a 


Ha:a^=0  to  that  of  the  corTespondmg  likalihood  ratio(LR)-test  for  different  degrees 
of  imbaliuire.  They  found  that  the  LR-tost  could  be  appreciably  more  powerful  than 
the  F-test  for  extremely  unbalanced  deeigns. 

Singh  (1991)  studied  the  effects  of  heteroscedaaticily  and  data  imbalance  on  the 
power  of  the  usual  F-test  of  WjroJwO  when  the  error  variances  in  the  one-way  random 
model  ere  heterogeneous.  He  found  that  beteroscedaaticity  increases  the  power  of  the 
tfflt  in  situations  where  the  larger  variances  are  associated  with  the  smaller  group 
sizes.  If  the  larger  variances  are  associated  with  the  larger  group  Kzes,  the  power  is 
decreased.  A small  increase  in  the  power,  doe  to  bcterosccdasticity,  was  noted  with 
balanced  data. 

El-Basiouni  and  Seely  (1096)  compared  demgns  on  the  basis  of  power  of  teste. 
Tliey  proposed  the  modiBed  harmonic  mean  test  procedure  for  variance  components 
and  compared  this  lest  with  the  locally  most  powerful  (LMP)  teat  and  with  the  Wald 
l«t  for  the  one-way  random  model.  They  considered  ratios  of  empirical  powers  of 
these  tests  to  the  power  of  the  most  powerful  (MP)  t«t  at  17=0.1, 0.3.  and  O.5.,  where 
17  = 0.  They  argued  that,  when  there  is  no  a priori  knowledge  of  17  (or  p),  or  when 
the  design  is  badly  unbalanced,  the  modified  harmoaic  mean  test  is  recommended. 
Since  powers  are  functions  of  p and  the  design  (namely,  k,  and  Af  for  Che  one-way 
dassiheation  model),  the  effects  of  demp  and  p on  the  powets  of  the  tests  should  be 
examined.  We  can  possibly  accomplish  this  by  establishing  an  empirical  relatioiLShip 
between  the  power  of  a test  and  IV.  and  p using  Khuri's  (1996)  approach. 

Thomas  and  Hultquisi  (197S)  studied  interval  estimation  of  vtkriance  components 
using  the  one-w^  random  model.  They  considered  a statistic  which  Is  a function  of 
the  harmonic  mean  of  the  m’a  and  of  the  sample  variance  of  the  treatment  means. 
They  argued  that  the  exact  distribution  of  this  statistic  can  be  approximated  by  a 
chi-squared  distribution.  The  adequacy  of  this  approximation  depends  on  the  desIp 


!iip  and  values  of  varlanc 


componsnis,  the  coverage  probability  of  a confidence  interval  for  variance  components 
was  obtained,  and  compared  with  the  nominal  value  of  the  confidence  coefficient,  l-o- 
Here  also,  we  can  investigate  bow  this  approximation  woriis  more  fnlly  using  Kburi’s 
(1996)  approach. 


CHAPTEH  3 
DESIGN  COMPARISONS 

3.1  Iptroduction 

The  quality  oI  estimation  of  variance  components  depends  to  a large  extent  on  the 
design  used  to  generate  the  response  date.  The  majority  of  the  published  work  to  this 
design  area  appeared  in  the  I9fi0's  and  1970's  [see  the  review  articica  Anderson 
(1975)  rtnd  Kliuri  and  Sahai  (1935)  for  a comprehensive  coverage  of  the  relevant 
literature].  For  the  most  part,  the  criteria  used  in  the  selection  of  dedgns  were  based 
on  single-valued  functions  such  as  the  trace  or  determinant  of  the  variance-covariance 

and  Gaylor,  1970), 

An  unsettling  feature  of  designs  for  vanance  components  estimation  is  their  de- 
pendenc,v  on  the  values  of  the  variance  components  themaelvea.  Thus  the  choice 
of  an  efficient  design  cannot  be  made  without  some  knowledge  of  these  parameters, 
specially  when  a single-valued  criterion  is  used.  The  same  applies  to  the  compari- 
son of  several  designs.  This  dependency  problem  was  noted  by  Senrle,  Caaella,  and 
McCulloch  (1992,  p.  76)  in  their  study  of  the  eifect  of  imbalance  on  the  variance  nf 
the  ANOVA  estimator  of  the  among-groups  variance  component  [or  an  cme-w^ 
raiidum  model.  Lohr  (1995)  remarked  that  the  dependency  of  a criterion  function  on 
the  unknown  parameters  makes  the  design  problem  a nonlinear  one.  To  resolve  this 
problem,  Anderaon  (1931),  Mukerjee  and  Huda  (1936),  and  Giovagnoli  and  Sebastiani 
(1939)  obtained  locally  optimal  designs  in  the  sense  of  Chernoffi  (1953).  This  approach 
requires  the  speciheatioo  of  some  prior  knowledge  of  the  variance  components.  As  a 

17 


18 

r€sult,  the  efftcieccy  of  the  gertereCed  deeigne  will  be  dependent  on  the  reliobillty  of 
the  posttilfttcd  prior  knowledge.  Another  approach  to  the  deeign  dependency  problem 
in  based  on  the  use  of  the  Bayeeian  methodology,  which  regards  variance  components 
an  random  variables.  The  latter  approach  requires  the  specification  of  a prior  dis- 
tribution for  the  variance  components.  Obviomsly,  different  prior  distributions  can 
result  In  different  designs.  Estimators  are  obtained  ^iii  posterior  distributions  of 
the  variance  components. 

Another  factor  that  affects  the  chedee  of  design  is  the  method  used  to  estimate 
the  model's  variance  components.  Mune  and  Anderson  (1978)  compared  designs, 
using  maximum  likelihood  (ML),  for  the  estimation  of  variance  comi>oiiems  in  a two- 
way  cross  classification  random  model.  Goldsmith  and  Gaylor  (1970)  used  analysis 
of  variance  (AXOVA)  estimation  to  compare  dealgns  for  a two-fold  nested  random 

In  this  chapter,  a graphical  approarh  is  presented  for  comparing  designs  for  esti- 
mating variance  comitoiienrs,  mt  well  a.s  for  com |>ariiig  different  methods  of  estimation 
with  a particular  design.  The  proposed  approarh  uses  the  so-ealJed  quantile  dispeiv 
slon  graphs  (QDGs)  that  were  first  Introduced  in  Khurl  (1097)  In  connection  with 
the  balanced  two-way  random  model.  The  QDGs  consist  of  plots  of  the  maxima 
and  minima,  over  some  region  in  the  parameter  space,  of  the  quantiles  of  a variance 
component  estimatur.  An  approach  using  the  so-called  empirical  quantile  dispersion 
graphs  (EQDGs)  is  proposed  In  this  study  when  exact  quantile  values  are  unavail- 
able. The  EQDGs  arc  obtaioed  by  plotting  the  maxima  end  minJmo  of  the  empirical 
quantiles  of  a variance  component  estimator.  These  plots  provide  a comprehensive 
picture  of  the  quality  of  estimation  with  a pertieular  design.  They  can  also  be  used  to 
compare  designs  for  the  seme  model.  Two  methods  of  ^timatiou  will  be  conridered, 
namely,  ANOVA  and  ML. 


[n  Section  3.2,  designs  for  one<way  random  model  are  considered.  Two-way  cross 
clateiiBcatioii  random  model  is  investigated  in  Section  3.3.  Section  3.4  is  devoted  for 
three-fold  nested  random  model,  and  is  followed  by  conclu^n  In  Section  3.3. 

3.2  The  One-Way  Random  Model 

3.2.1  Introduction 

Consider  the  imbalanced  one-way  random  model, 

(3.11 

I = 1,2,'  • • ,Ar;  y = 1,2,--  - ,ttn  where  /r  is  a fixed  unknown  parameter,  o,  is  the 
effect  of  level  i of  the  treatment  factor,  and  c,,  is  a random  error.  We  assume  that 
a,  and  ctj  ore  independently  distributed  such  that  a,  N[0,a^),  tif  — fV(0,o^), 
I = 1,2,--  - , and  j = 1,2,--  - ,n,.  The  model  Is  balnnccd  if  the  tic's  are  equal,  that 

is,  rii  = = - ' • = n*  = n.  Model  (3.1)  can  be  written  in  vector  fonn  as 

V = (i1a  -I-  (©^,1„,1 0-1-  e,  (3.2) 

where  N = Ef.i"!.  V = (llii.Pie,  • • • .Bi.,.- ••  ,l»i,»«3.- • • .Ift.,)'.  In  ‘he  vector 
of  ones  of  order  A^xl,  © is  the  symbol  of  direct  sum  of  matrices,  o = (oi,  oj,  - • • ,as)', 
and  6 is  defined  analogously  to  y.  The  mean  of  p is  pi  w,  and  the  variance-cuvariance 

E = <rj  ©f„  J,,  -hoj/w,  (3.3) 

where  fw  is  the  identity  matrix  of  order  N xN  and  J„,  is  the  mntrix  of  ones  of  order 

The  proposed  graphical  technique  for  comparing  designs  for  model  (3,1)  uses 
quantiles  of  the  distribution  of  where  denotes  an  estimator  of  o^.  ANOVA 
estimation  is  used  to  Sections  3,2.2  and  3.2.3  while  ML  estimatioD  is  considered  in 
Section  3.2.4.  Comparisons  of  these  cstimstlon  methods  for  given  designs  are  made 


in  Section  3.2.5.  In  Section  3.2.4,  deigns  for  estimating  the  ratio  p — gr??  ^ 
compared  usng  also  ANOVA  and  ML  estimation. 

3.2.2  Quantiles  of  ^ mnmtANQVA  eetitnation 
The  ANOVA  table  for  model  (3.1)  ia  of  the  fonn 

Source  d-f.  Sum  of  Squares  Expected  Mean  Squares 
Among  groups  1;  - 1 SSa 

Within  groups  N — k SSs 

Total(adjusled)  N-l  SSr 


S5„  = v'Qv.  (3.5) 

SSb  = y'Ry.  (3.6) 

Q = ffi‘„  - ij„,  (3.7) 

R = (3.8) 

The  matrix  QS  can  be  expressed  as  (sec  Khuri,  Mathew,  and  Sinha,  1998,  p.  54) 

QS  = ^©?.i  Jn,  - ^Ixd'J  aj  + I®?.,  i Jx]  of,  (3.9) 


where  d = (nil„,  ; : rifcl„j ).  In  general,  this  matrix  is  not  a multiple  of 

an  idempotent  matrix  unless  the  design  is  balanced.  Hence,  55a  is  not  distributed  as 


» sc&l&r  miiliiple  of  a chi-squared  variate  as  is  the  ease  with  balsnred  models  (see,  tor 
example,  Seerle,  Casella,  aud  McCulloch,  1902,  Theorem  S,2,  Appendix  S,5,  p,  467], 
The  tollowing  theorem  shows  that  SSa  can  be  expressed  as  a linear  combinatioD  of 
independent  central  chi-squaretl  varintes. 

Theorem  3.2.1  For  n p-dimenaional  vector  x ~ ,V(0,  S)  and  asymmetric  matrix  A. 
X Ax  can  be  expressed  as  a linear  combination  of  independent  central  chi-squared 
variates  of  the  form 

where  A|,Av,  - - ,A,  are  the  distinci  nonzero  eigenvalues  of  AD  with  multiplicities 
' ■'  > Vr,  respectlvTly- 

Proof  (See  Johnson  and  Kotz,  1970,  Ch,  29,  pp.  150-153).  Q 

Corollary  3.2.1  The  result  of  Theorem  3.2.1  remains  true  if  x ~ jV(si,S)  provided 
that  ;t' Ap  - 0 and  A is  symmetric  positive  semidehnite. 

Proof  .Assuming  x --  .V(*t,E),  let  y = E"'x  - iV(E"J/i,/,v).  Then,  by  the 


xAz  = p'siAZiy 
= yPAP'y. 

where  A - diag{6,),  i s 1,2,  > > > , (V,  with  J,  being  the  eigenvalue  of  A£^,  and 
P is  an  orthogonal  matrix  whose  columns  are  ortbonormal  eigenvectors  of  S^AS^ 
corresponding  to  the  J,’s.  Let  ta  = P'y,  then  to  N(P'S~^ti,2)  and 


e.  = where  P.  is  a matrix  of  order  W x v.  whose  columas 

are  eiRenvectore  corrrapondiog  to  A,.  Thus,  x Az  & E^.A,  if  sad  only  it 
ft  = 0 for  all  1 = 1,2, --.r,  lliat  is,  ^'J:-!P,P;E-Sk  = 0 for  all  i = 1,2,  • • • ,r,  or 
eqaivaleatly,  if  and  only  if  Jt'E-*(E^i  = 0 (the  sufficiency  is  due  to 

A,  > 0 which  ia  true  if  A is  positive  semidefinite),  that  is,  ,*  E-IPAP'E-i>i  = D, 


(3-10) 


and  Q is  the  matrix  defined  in  formula  (3.7). 


Proor  The  expected  \ulues  of  the  mean  squareti  of  the  between  group  effect  and 
the  error  tenn  are  £(MSo)  = d aj  + o?  and  £(AfSa)  = respectively,  whore 
j s I^[iV  - j;  n?).  Hence,  the  ANOVA  estimator  of  <rj  is 

nix  = 5l.WS.-Wsl 
3S.  SSe 

Note  that  Q is  a symmetric  positive  semidriinite  matrix.  Since  ^(v}  s ply  and 
Qly  = 0,  by  Corollary  3.2.1, 5S.  = y'Qv  - HJ.i  AjJtS,,  where  A|,Aj,  ••  • ,A.  are  the 
distinct  positive  eigenvalues  of  QE  with  multiplicities  mi.mj,--'  , m,,  respectiv.ely, 
and  E is  defmed  as  in  equation  (3.3).  Let  E'  = ^ Then 

(311) 


where  AJ,  Aj,  • • • ,A;  are  the  distinct  positive  eigenvalues  ofQEV  On  the  other  hand, 
the  error  sum  of  squares  is  distributed  as 

SSE~ahl-k.  (3-12) 

independently  from  5S..  By  combining  equations  (3.11)  and  (3.12),  we  find  that 
is  distributed  as 


Cnrollarv  3.2.2  For  the  balanced  one-way  random  model,  y,j  = p -H  ' 
1.2,  j = 1,2,  ••  ,n,  the  scaled  ANOVA  estimator  of  oj.  W’„,,  = 
trihuted  as 


(fc-l)nq 


(3.14) 


Proor  In  Ibis  c»se,  the  expected  values  of  mcar  squares  of  between  group  effect 
nnd  the  error  effect  are  E{MS.)  = na‘  + o’,  ami  E{MSe)  = <r?.  Therefore,  the 
ANOVA  estimator  of  oj  is  pven  by  oj^  = ;(AfS„  - .WS«).  Since  for  balanced 
model,  SSo  ~ (noj  + oflxf.n  and  d in  (3.4)  is  n, 

-S  elmr; + ■>,')  » _ _f? j (3.15) 


Hence,  the  result  foUowe. 


CoroUaitiSJ  The  probability  of  a negative  ANOVA  estimate  of  oj  for  the  unbal- 
anced one-wav  random  model  is 


P{Sij 


' rbE*: 


For  the  balanced  ono-way  random  nmdel,  it  la 

Proof  These  formulas  follow  from  (3.10)  and  (3,14),  respectively.  D 

The  exact  distribution  of  Wai  in  (3.10)  and  (3.14)  can  be  determined  by  its 
quantiles.  The  />“  quantile  of  iV)^,  denoted  by  g{p),  is  such  that  PlH',a  < q(p))  = 
p.  Since  li',x  is  a linear  combination  of  independent  chi-squared  variates,  an  exact 
tabulation  of  its  cuitiulative  distribution  function,  F[WaA  < d(p)l-  can  be  easily 
obtained  usbg  Davies'  (1980)  algorithm.  Tbe  algorithm  itself  (identified  as  AS155) 
is  accessblc  via  STATUE,  an  E-mail  and  file  transfer  protncol-based  retrieval  system 
tor  statlstica)  software.  A tabulation  of  the  quantiles  j(p)  for  sJven  values  of  p is 
therefore  posable  through  Davies'  algorithm. 


The  if'  quintile,  in  (3.10)  depends  on  the  unknown  value  of  ij  = ^ 

as  well  as  oo  the  design  used  by  virtue  of  (3.0).  By  a design  we  mean  a speciftcstion 
of  the  values  of  R|,ne,  * " .n^  for  a given  h in  the  model  of  (3.1).  Such  a de^gxi  is 
denoted  by  D.  Thus,  D = (ni.na,.--  ,nj}.  The  dependency  of  thep“  quantile  on  i) 
and  D will  henceforth  be  indicated  by  writing  9o(p,7). 

Khuri  (1997)  intioduced  the  notion  of  quantile  dispeision  graphs  (QDGs)  as  a 
means  to  compare  designs  for  estimating  variance  components  on  the  basis  of  tbeir 
quantile  profilts.  For  model  (3.1),  the  QElGs  condst  of  plots  of  Qg“(p)  and  Q5“(P) 
against  p.  where 

Qo“(>’)  = (3.18) 

Co^(p)  = nim8o(P.'!).  (8.19) 

Here,  denotes  a region  in  the  parameter  space  of  q.  Such  plots  provide  an  overall 
pictorial  assessment  of  the  quality  of  A140VA  estimation  for  a pveii  design  D.  Flat 
plots  and  low  values  ofQg“(p)  -<Jo""(p)  l8e  range  of  p Indicate  a high  quality 

of  estimation  and  are  therefore  desirable.  This  implies  stability  in  the  values  of  the 
AN'OVA  estimator  of  oj,  and  little  sensitivity  to  changes  In  the  v^ues  of  q.  In  other 
words,  one  design  is  preferred  over  another  if  the  QDGs  of  the  former  display  smaller 
changes  in  the  quantile  values  than  those  for  the  latter. 

Let  us  consider  the  followtug  two  designs  for  model  (3.1)  with  k = 5 and  N = 25: 
Du  = {1,1,1,11,11}  (3.20) 

D^  = {1,1,1,1,21}.  (3.21) 

These  designs  were  initially  conddered  by  Searle,  Casella,  and  McCulloch  (1992,  pp. 
76-76}  in  their  study  of  the  effect  of  imbalance  on  the  variance  of  Oas*  ^ addition  to 


Dgi  and  lat  ua  also  consider  the  balanced  design  for  ^ - b and  N ~ 25,  namely 


Dj  = {5.5, 5.5, 5).  (3-22) 

Tbe  estimation  capabilities  of  these  three  designs,  with  regard  to  9^,  will  now  be 
compared  using  the  corresponding  QDGs  of  in  (3-10)  and  (3,14), 

For  the  balanced  design  by  (3.14) 

W 


Applying  now  (3.10)  in  conjunction  with  D^,  we  get 

Wi*i +*is>ti +*is*l)  (3-24) 

where  Ajj,  A),,  and  Ajj  are  the  distinct  nonaero  eigenvalues  of  in  Tlieorcm  3.2.2 
with  multiplicities  1,  1,  and  2,  respectively.  Similarly,  for  the  design  Z>us,  we  have 
»;4S±(A5,A?  + Ai;d)-3^Ai„  (3.25) 

where,  in  this  case,  Aj,  and  Aj,  arc  the  district  nonzero  eigenvalues  of  QE'  with 
multiplicilies  1 and  3,  respectively. 

The  quantiles  of  IFoa  for  designs  Di,  D^ir  and  0^2  are  now  computed  using 
formulas  (3,23),  (3.24),  and  (3.25),  respectively,  iu  conjunction  with  Davies'  algorithm 
as  was  mentioned  enrlier.  For  each  design,  values  of  and  Q^'^'fp)  in  (3.18)  and 

(3.19)  are  obtained  for  each  of  several  values  of  p (0  <p  < 1).  Here,  the  optimization 
of  en(p,h)  witli  respect  to  i;  is  carried  out  over  a region  S„  of  the  parameter  space 
given  by 


S,  = (i?  I 0.2  S 1 < 10.0). 


(3.26) 


More  apedficaJIy.  for  a given  p,  qoij3,v)  computed  using  the  following  values 
of  IJ  : I)  = 0.2(0.1)10.0.  The  maidinuni  and  minimum  of  the  resulting  values  of 
«o(p.  o)  provide  approximate  values  of  (3o“Cp)  and  Q^(p).  The  results  are  given 
in  Table  3.1. 

The  QDGa  for  Dt,  H.i,  and  are  conatrueted  on  the  basis  of  Table  3.1.  A 
display  of  the  QDGs  lor  each  pair  cf  designs  is  given  in  Figure  3.1.  The  following 
observations  can  bo  made: 

(a)  The  variability  in  the  quantile  values  with  respect  to  q for  £>i  and  £l.i  are  quite 

similar  [see  Figure  3.1(a)].  Design  D^,  however,  displays  a larger  degree  of 
variabiUty  chan  both  Dt  and  D,i  [see  Figure  3.1(hx:)].  This  holds  true  for  all 
values  of  p. 

(b)  There  is  a greater  degree  of  uniformity  in  the  quantile  values  for  Dt  and 
than  for  0„3.  Thus  more  stadlity  in  the  values  of  is  achieved  with  Dh  or 
D,i  than  with  D^. 

On  the  basis  of  these  observations  it  can  be  easily  concluded  that  the  estimation 
capability  of  Dbi  is  better  than  that  of  D^.  The  performance  of  Dy\  is  surpridngly 
similar  Co  that  of  Di  with  the  latter  being  slightly  better  than  the  former  [see  Fig- 
ure 3.1(a)).  Thus  the  use  of  the  QDGs  made  it  quite  easy  to  discriminate  among  the 
three  designs,  particularly  between  D„i  and  D^.  Such  a task,  however,  is  not  always 
possible  with  other  approadies.  For  example,  using  the  value  of  as  a crite- 

rion for  comparing  demgns,  which  Searle,  Casella,  and  McCulloch  (1992,  pp.  75-78) 
used,  it  is  not  possible  to  decide  which  of  D^i  and  Ao  is  the  better  design  uniformly 
over  all  values  of  q.  as  can  be  seen  from  Table  3.2.  They  noted  that  for  some  values 
of  q,  Ai  is  better  than  Dns  becauae  A,  has  smaller  values  of  than  D^s, 
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(a)  DJ)  vs  Djji  (b)  OJj « Dji2  (c)  0_ul  vs  Oji2 


Figure  3.1.  QDGs  of  Woa  fur  the  three  designs,  Dt,  D^\,  ■ 


Tible  3.3.  Values  of 


lor  designs  D>,  D,i,  and 


D,  = {5,5, 5, 5, 5)  0.1053 

= 0.1765 


0.7240  2.4240 
1.4114  4.9945 


D.s  = (l,  1,1, 1,21} 


0.6620  1.5370  4.4815 


10 

52.0240 

113.3360 

85.3704 


However,  focoiher  values  of  i|  (i)  > 1),  Dgi  appenretobc  worse  than  D^.  They  indi- 
eatod  that  “The  value  of  i}  at  which  Var(^)  changes  from  increaang  to  decreasing 
as  unbalancedneas  increases  is,  of  course,  not  the  ewme  for  all  situations,  but  depends 
in  no  simple  way  on  JV,  k,  and  the  n,<vaJuea"  (Searlc,  Casella,  and  McCulloch,  1992, 
p.  76).  However,  based  on  the  Figure  3.1,  we  can  easly  see  the  preference  of  design 
fl„i  to  demgn  D^,  and  the  similarity  between  designs  Dt  and  fOui. 

.1.2.4  Comi>arison  of  desisns  usine  ML  estimation 


In  this  section,  we  reconrider  the  problem  of  detngn  comparison  for  model  (3.1) 
u»ng  maximum  likeliliood  (ML)  estimation  of  o^.  Let  denote  the  ML  estimate  of 
and  let  tV'^  = Unlike  the  case  of  ANOVA  estimation  in  Section  3.2.3,  6^^, 
does  not  have  a closed-form  expreslon  when  the  dataset  is  unbalanced.  The  estimate 
°tily  be  computed  numerically  using,  for  example,  PROC  MIXED  procedure 
in  SAS  (1997),  Furthermore,  the  distribution  of  is  known  only  asymptotically. 
Ctmsequently,  the  qurmtiles,  and  hence  the  QDGs  for  iVou,  cannot  be  olttained  ex* 
actly.  Thtty  can,  however,  be  obtained  approximately  using  the  so-called  empirical 
(|iiantile  dispersion  graphs  (EQDGs)  for  These  can  be  developed  througli  Monte 
Carlo  sirmilation.  An  outline  of  this  process  follows: 
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(«)  Without  liny  loss  of  generality  we  Cin  assume  that  ji  = 0 in  model  (3.1).  We 

il.,  = £.,  + <«,  . = l,2,---,k;  j = l,2,  -,n..  (3.27) 

We  recall  that  o — fh*  Hence,  a,  W(0,ii^j)  and  Cy  iV(O.of).  Let  U denote 
a rc^on  in  the  parameter  space  of  the  pair  {rj,  of  )• 

(b)  For  a given  (t|,(rj)  in  U and  a given  design  25  = (ni.nj,  • • • ,n»},  a random 
vector  y is  generated  on  the  basis  of  model  (3.27):  one  random  number,  which 
we  denote  by  t,,  Is  generated  from  W(0,  osf ).  In  addition,  n,  independent  ran- 
dom numbers,  denoted  ty  Sn.Sn.- **  .Sin.,  are  generated  from  ^(0,rrf)  inde- 
pendently oft*.  This  process  produces  the  vector  (ti-hsii,ti-i-sij,  ■-  - , Ir  + Sin,)', 
which  represents  a sample  of  response  values  from  Che  i*^  group  (t  s 1,2,  - • • ,k). 
By  repeating  the  procets  independently  k times  and  combining  the  results,  we 
obtain  a vector  y consiating  of  n,  = hf  observations. 

(c)  Lsing  tbe  response  vector  y obtmned  in  step  (b),  a ML  estimate,  of 

is  computed  using,  for  example,  PROC  MIXED  procedure  In  SA5  (1297).  A 
value  of  W'o*r  — ^ is  then  obtained. 

(d)  By  repeating  steps  (b)  and  (c)  a sulficlent  number  of  times,  we  can  obtain 
empirical  quantiles  of  CVair  for  selected  values  of  p for  a given  (q.ec).  We 
denote  the  empirical  quantile  thus  obtained  by  qo(p,q,o^), 

(e)  Steps  (b),  (c),  and  (d)  are  repeated  several  times  corresponding  to  different  values 

of  (q.of)  selected  from  the  re^un  V.  The  set  of  such  selected  pairs  is  denoted 
by  T.  Thus  T C If. 


of  the  cnipirical  quantile  values  over  the  set  T 


in  step  (e)  arc  obtained.  We  then  have 

^3“(p)  = (3.2S) 

(jrw  = (3.29) 

(g)  Vaiucs  of  Qo"(P)  “<1  Qo"'(p)  mc  plotted  against  p to  obtain  the  EQDGs  for 

IV.«. 

Let  us  now  illustrate  the  application  of  the  above  procedure  uung  the  same  de- 
dgna,  Dh,  D^u  and  as  in  Section  3.2.3.  The  region  V in  chosen  as 

= {(h.  ir?)  I 0.2  < q < 10.0,  0.2  < of  < 2.0} . (3J0) 

Tiie  set  T ineniioiied  in  step  (e)  condsts  of  300  points  selected  from  U such  Uiat 

T = {(O.ff?)  I 0 = 0-2(0J)10.0,  of  = 0.2(0.2)2.0>. 

A SAS  Macro  (SAS,  1990)  wan  written  to  compute  9o(p.0.i^)-  A total  of  2000 
replicatione  in  step  (c)  were  made  to  calculate  values  of  dpfp.q.of)  for  each  design 
and  each  (q.of)  in  the  region  T.  Thbic  3.3  gives  values  of  Qp“(p)  and  (3o^''(p)  In 
(3.20)  and  (3.29),  respectively,  forila,  and  ^^3.  Plots  of  these  values  for  each  pair 
of  designs  are  given  in  Figure  3.2.  As  in  Figure  3.1,  we  note  hero  that  the  variability 
in  tbe  quantile  values  for  Dt  and  £>„i  ore  similar  (see  Figure  3.2(a)],  However,  the 
superiority  of  D,i  to  D„j  is  not  as  evident  as  in  Figure  3.1  (see  Figure  3.2(c)|. 

3.2,5  Copiaarison  of  AMOVA  and  ML  estimatioo 

We  can  use  the  information  obtained  in  Sections  3.2.3  and  3.2.4  to  compare 
ANOVA  estimation  with  that  of  ML  with  regard  to  irj  in  0 limited  sense  sbice  tbe 
reruns  S,  and  V in  (3.26)  and  (3.30),  respeetively,  ere  not  the  same.  The  QDGs 


SSS.S  sr=rc'“ 


(a)0_6ve0jl  (b)D_bv5  0jj2  (s)  D_u1  « D_u2 


Fifiiri!  3.2.  EQDGb  of  IV.m  for  Uie  three  deigns,  Di,  D.i,  and  D, 


for  If'oj  and  the  EQDGs  for  H'.u  are  shown  together  for  each  design  in  Figure  3.3. 
For  the  balanced  design  D^,  the  graphs  for  ANOVA  and  ML  are  almost  identieal, 
except  for  the  pnsaiWIity  of  getting  negative  values  with  the  ANOVA  estimate  |sse 
Figure  3.3(a)].  Note  that  for  iudanced  daU,  restiicied  maximum  likelihood  (REML) 
estimates  are  identical  to  the  ANOVA  estimates  when  the  latter  arc  nonnegative. 
For  design  D„,,  the  similarity  between  ANOVA  and  ML  values  is  moderate  [see  Fig- 
ure 3.3(b)).  The  frequency  of  getting  negative  ANOVA  estimates  is  greater  udth  D^\ 
than  with  Finally,  for  design  the  distinction  between  ANOVA  and  ML  Is 
more  pronounced  than  for  the  other  two  deigns  [see  Figure  3.3(c)].  In  this  case,  ML 
is  preferred  over  ANOVA  in  terms  of  etobility  in  the  quantile  values. 

3 2 6 Estimation  of  the  inlracliss  correlation  coefficient 


The  intraclass  correlation  coefheient  (ICC)  for  model  (3.1)  is  defined  as  the  ratio 
p ~ . This  is  a measure  of  association  among  the  responses  from  the  same 

treatment  group.  We  now  discuse  the  cempariaon  of  demgns  for  estimating  p. 

The  ANOVA  estimate  of  p la  pa  = where  dja  ^<a  *te  the  ANOVA 

estimates  of  and  a},  respectively.  Similarly,  the  ML  estimate  of  p is  Pvr  = , 

where  Slu  *?«  the  ML  estimates  of  oj  and  of,  respectively.  The  quantiles 
of  pa  can  be  obtained  by  using  Davies'  (1960)  algorltbm;  If  qf{p)  denotes  the 
quantile  of  pa.  then 

P[ba  < »a(p)]  = P|{1  - 9,(P)}^  -9,(p]^?a  < 0].  (3.31) 


Since  ofa  - ^ and  dja  P'®'  *>y  (3-13),  we  have  on  the  right-hand  ride  of  (3.31) 

applies.  On  the  other  hand,  the  quantiles  of  par  can  be  obtained  by  rimulation  as  In 
Section  3.2.4. 
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Figure  3.3.  QDGs  of  (i'o.4  find  EQDGs  of  (i’oAr  for  Che  three  designs.  Dt,  D^i,  and 

D.1 


It)  0_0  »s  DJ|1  (b)  DJ  vs  Djj2  (e)  D_ul  vs  OjiS 


Fiipire  3.4.  QDGs  of  Pa  for  the  three  deigns,  Dt,  and 


(U0J)«Djj1  ffl)DJivsO_ll2  (a)D_u1«D_u2 


Figure  3.S.  EQDGe  of  pu  for  i 


uid 


rtspcclively.  Wc  note  that  the  two  estimation  oietbods  do  not  exhibit  any  remnrk- 
nblc  diflerences  in  cbe  pbts,  cxrcpt  for  the  pontibilicy  of  haviiii:  negative  ANOVA 
estimates,  especially  with  design  0^2- 


(a)  0_b  lei  D_iii  (e) 


Figure  3.6.  CJDGs  of  flj  and  EQDCb  ot^M  for  the  three  derngna,  Di,  D,i,  and  O.j 


3.3  The  Two-Way  Cro86  qaaification  Random  Model 


3.3.1  latrpductioD 

Consider  the  unbalanced  two-way  cross  classificalion  without  interaction  random 
model, 

= (3.32) 

I = 1,2.  • • • .r;  j = 1,2,  • • • , s;  k — 0, 1,2,--  - ,n,^,  where  tr  is  a fixed  unknown 
parameter,  o,  is  the  effect  of  the  level  of  the  row  factor  A,  ^ is  the  effect  of  the 
level  of  the  column  factor  B,  and  Cys  <e  a random  experimental  error.  We  asuine  that 
a,  and  3,  are  independently  distributed  such  that  a,  ~ A'(0,o^)  and  B,  ~ fV(0,oJ), 
I = 1,2,  • • • , r,  j = 1,2,--'  ,a,  respectively.  Moreover,  eys  is  distributed  as  W(0,  of), 
independently  from  a,  and  k = 0, 1, 2,  • • • ,n,^.  The  model  is  balanced  if  the  ny’s 
are  equal,  that  is,  nu  s n^  s • • • s 0,1  - n (,£  0).  Model  (3.32)  can  be  writteo  in 
vector  form  as 

V = f.l«-tX,e.-(-X,fl-t-e,  (3.33) 

l,v  = the  vector  of  ont«  of  order  IV  X 1, 

-2fi  — where  n,.  = ^ riy, 

0 = (/5t,y3s,---,fl.)’, 
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and  € ia  defined  anologoiiflly  to  tr<  an  /V  x 1 vector  of  oheervncioiut.  The  mean  of  v ia 
and  ita  vaHance<cuvariaDce  matrix,  denoted  by  S,  is  oF  the  Form 

S = XiX',al  + XtX',a\  + (3.34) 

Farther,  let  us  define  the  following  matrices  For  l^er  use. 


Z,  = ll„  ; X,  : X,]. 

(3.35) 

P,  = z,(z',z,)‘zi.  and 

P,  = Z,(Z',Zj)'Zi, 

where  A~  is  a generalized  Inverse  oF  a matrix  A.  Note  that  Pq  = s-fw,  where 

Fn  this  section,  we  conrider  the  quantiles  oF  the  distributions  oF  ond  ^ for 
comparing  tiesigns  as  well  as  comparing  two  estimation  methods  of  the  variance  com* 
ponents,  namely,  ANOVA  and  ML,  where  and  denote  estimatora  of  ^nd  o|, 
respectively.  Three  designs  are  cliosen  for  comparison  purposes.  These  designs,  which 
were  used  by  Muse  and  Anderson  (197S),  namely,  B*3d,  Bf32>36,  and  Of33‘35,  will  be 
deRned  in  detail  in  Section  3.3.3.  ANOVA  estimators  are  investigated  in  Sections  3.3.2 
and  3.3.3,  while  ML  estimators  are  considered  in  Section  3.3.4.  Comparisons  of  these 
estimation  methods  for  eadt  design  arc  mode  in  Section  3.3.5. 

3.3.2  QuaciUeaof^  and  ^ uaiM  ANOVA  estimation 

We  shall  use  Type  f sums  of  squares  (SS)  to  measure  the  contribution  of  factors 
A and  B.  Type  I SS  ore  sometimes  referred  to  as  sequential  sums  of  squares. 


An  ANOVA  table  For  model  (3.32)  is 


Source 

d.f. 

Type  I Sum  of  Sqtiares 

Expected  Mean  Squares 

A 

/i 

SSa=fi(o|p)=fl(p,Q)-«(p) 

fc|tr^  + kiOg  -h  of 

B 

/j 

SSg=R(fl|p,  i»)=fl(p,  tt,  fl)-R(p,  q) 

*t»s  + of 

Error 

I, 

SSc=y’v-R(p,“i3) 

of 

Tht.(adj.) 

N-l 

SSt 

where  /i,  /j,  aod  &re  the  numbers  of  degrees  of  freedom  for  A,  B,  end  ihe  error 
term,  reeprccfvely.  More  epecihcaily,  fi  and  fj  are  the  numbers  of  linearly  indepen* 
dent  Type  I eetimable  functions  of  the  factors  A and  B,  respectively.  A function  of 
parameters  Is  defined  as  estimable  if  there  is  a linear  function  of  the  vector  of  obser- 
vations y whose  expected  value  is  identical  to  tits  function  of  the  parameters.  Type 
I esrimabte  function  is  an  estimable  Function  based  on  the  partitioning  of  Type  J S5. 
For  a more  rigorous  deSoitlon  of  ao  estimable  function,  see  Graybill  (1076,  pp.  4S4- 
486).  Also,  /s  = A'  - — /j  - 1.  The  notation  fi(*)  denotes  e reduction  io  the  sum 

of  squares  with  the  contents  of  the  brechets  indicating  the  model  fitted  (see  Senrle, 
1971,  pp.  246-249).  For  example,  il(p,  a,  d)  is  the  reductioo  in  sum  of  squares  for 
fitting  model  (3.S2).  Similarly,  R{ii,  a)  is  the  reduction  in  sum  of  squares  for  fitting 
the  model  = p 4*  a,  -h  Sgs-  Hence  R{/A.a,d)  — A[p,a)  ta  the  reduction  m sum 
of  squares  due  to  fitting  a model  coutaining  p,  an  o-factor,  and  a if-Cactor  having 
already  titled  one  with  p and  an  o-factor.  This  Is  a measure  of  the  extent  to  which 
a model  can  explain  more  of  the  variation  in  y,js  by  including  a column  effect  (d)  in 
addition  to  an  overall  effect  (p)  and  a row  effect  (a).  The  cpefficieota  of  the  variance 
components  in  its  expected  mean  squares,  namely,  kj , k^,  and  hj,  can  be  obtained  by 
using,  for  example,  the  synthesis  method  of  Hartley  and  Rno  (see  Hartley.  1967). 


I 55xi  SSst  &nd  SS^ 


SSa  = yQ^y.  (3.36) 

SSb  = vQ^V,  (3.37) 

SSb  = yRy,  (3.38) 

Qa  = Pi  - Pa,  (3.39) 

Qb  = Pj-Pi.  (3.40) 

H = (3.41) 


uid  the  matrices  Pg.  P„  and  P,  are  defined  io  (3.35).  The  mean  squares  are 
denot4?d  by  AfS<,  A/5e,and  MSe,  respectively.  The  sums  of  squares,  SSa  andSSs, 

model,  since,  with  S in  (3.34),  the  matrices  and  Q^S  are  not,  la  general,  a 

Theorem  3.2.1  and  Coiollary  3.2.1  Co  express  SSa  and  SSg  as  linear  combinations  of 
independent  central  ohi^squared  variates. 

Let  us  denote  Type  I ANOVA  estimators  of  and  by  o^a  *tid  respec- 
tively. Let  us  also  denote  the  scaled  vereion  of  these  estimatoca  by  W„a  = ^ and 
WiA  = respectively.  The  distributions  of  W,a  and  W,a  can  then  be  obtained 

Theorem  3.3.1  For  the  unbalanced  two-way  cross  clandcaCion  without  interaction 
random  model  (3.32),  the  scaled  Type  I ANOVA  estimatoca  of  o‘  and  o|.  IF,.,,  and 


Wit,  nspcnively,  are  distribuled 


fi  E-rjX,’  Slid  {3.42) 

S (3-43) 


g — number  of  diatinct  po&tive  eigenvalues  oF  Qa£n, 
h - number  of  distinct  pcsitive  eigenvalues  of  QgSe, 

1i  - j“  positive  eigenvaiuesof  with  multiplicity  iy,  j = 1,2,-.  • ,9, 
6]  = pomtive  eigenvalues  of  QjSs  with  multiplicity  mj,  j = 1,2,  •••  ,h, 


Q, 

Oa 

E. 

Es 


*1 


= Xinr;  + x^x,  + 1„  . 


(3.44) 

(3.45) 
(3.4C) 
(3.47) 


The  matrices  Qt,  Qs,  and  R are  deBned  in  (3.39),  (3.40),  and  (3.41),  respectively, 
and  /i,  /e,  /g,  iri,  k-j,  and  kg  are  defined  in  an  ANOVA  table  which  is  given  earlier  in 
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si.  . 


crolUo-  3.2.1,  ~ n..3,;d,.  ''hero  Sui„-.S>.  are  U«  distinct  pooiUv 

rd  aa  ta  (3.34).  By  letting  Sj  = = X,X\  {^)  + X,X‘,  + In  {^),  th 


dlatrlbutl™  of  = obtained. 
The  Type  1 ANOVA  estimator  of  <r»  is 


E.  = = X,X\  * X,X,  (g)  + 1„  (g),  the  dii 


of  ^ is 
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Corollarv  3.3.1  For  Che  hatHCired  two-way  cross  clasUicatlon  without  ioCcractiori  ran- 
dom model,  yij  — /j  -e  tw  + + et;,  i — j = 1,2,  H'o,^  — and 

lies  = ^ an  distributed  as 


[r-l)ath*"'  (r- l)(e- l)»i)i 


Was  = 


Ly> ‘ V» 

(r-lKe-Drrb*"-"'-" 


(3.48) 

(3.49) 


Proof  In  this  case,  £(M5s)  = soj  + of,  £(MSs)  = roj  + of,  and  E{MSb)  = of. 
Since,  for  the  balanced  model,  S5s  (sof  -h  of)  ;l^i,  SSg  (roa  -h  of)  jf-j. 
5Sa  of  *4td  55s,  55s,  and  SSs  are  iodependetitly  distributed,  the 

ANOVA  estimators  of  of  and  of  are 


ofs  = -[M5s  - AfSal 
o|s  = ilMSs  - A/5aj 


sof  -h  of  , 
rai  + af  , 


(r-l)(.-i).*'"'>'-‘'' 

(r-I)(o-l)r*ir->M.-.)' 


resitectively.  Hence  the  result  follows.  Q 

The  exact  distribution  of  IV.s  in  (3.42)  and  (3.48),  and  W,,  in  (3.43)  and  (3.49) 
can  now  be  determined  by  their  quontilee.  We  con  use  Davies'  (1980)  algocichm  to 
obtain  exact  values  of  quantiles. 


The  p“  quantile  of  ll'os,  denoted  by  q(p),  is  such  that  P{WaA  < 9(p)l  = P-  For 
r model  (3.32),  this  quantile  depends  on  the  unknown  values  of  (h  = nndih  = ^ 


well  &s  cm  the  design  used  by  virtue  of  (3.44)  end  (3.4S).  The  secie  epplies  for  the 
if'  (luatitile  of  K'm  with  (3.45)  and  (3.47).  By  a design  we  mean  a spedfleation  of 
the  values  of  tzn,  H|3,  • . . .rv,  for  a given  r end  e in  the  model  of  (3.32).  We  denote  a 
design  ty  D.  The  dependency  of  the  jf  quantile  on  pi , Ps,  and  D will  therefore  be 
indicated  by  writing  W„a.  and  ’ll, ’}>)  for  Wee- 

Let  ufi  denote  a region  in  the  parameter  space  of  ttl  and  by  5(3,^).  The  QDGs 
for  W„A  are  obtained  by  plotting  0S2(p)  and  OgStP)  “gainst  p.  where 


0S?(P) 


(3.5D) 


(3.51) 


Similaify,  the  QDOs  for  are  obtained  by  plotting  Qo^(p)  end  Q"e(P)  “gainst 


Muse  and  Anderson  (1973)  compared  several  designs  for  model  (3.32).  Among  those 
designs,  the  following  three  designs  are  of  particular  intertst: 

(1)  the  balanced  design  with  r s 5,  5 a 6.  and  W w 36,  which  we  denote  by  L>is; 

(2)  the  so-called  hofonced  dr^oinl  nxtangl^s  design  with  r w 13,  s = 18,  and  N = 36, 

which  we  denote  by  Duj.  The  incidence  matrix  for  this  design  is  X = ItSJt. 
where  9 is  the  symbol  of  Kronecltcr  product  of  matrices,  namely  B = 
(Ui]£},  where  o,;  is  the  (t,  j)“  element  of  the  matrix  A. 

By  an  hicidenee  matrix  X,  we  mean  the  r x a matrix  having  elements  that  are 
all  0 or  1.  The  location  of  the  O’s  and  I’s  represents  the  incidence  of  terms  of 


the  model  correspoodiog  to  the  celli,  i = 1,2,  • • • ,r,  j = 1,2,  • • • , s,  of 
the  cUte. 

(3)  the  BO-ealJed  off-dta^oTiai  design  with  r = 18,  s = 18,  and  A'  = 38,  which  we 
denote  by  D^.  The  incidence  matrix  for  this  deagn  is  JC  = /«  » (Jj  - /j). 
Muse  and  Anderson  (1978)  compared  these  designs  using  asymptotic  variance  - co- 
variance  matrix  of  ML  estimators  of  variance  components.  For  example,  they  used 
the  trace  of  the  matrix,  Var(^),  and  Var{6}),  where  dj  and  dj  are  ML  estimators 
of  rr^  and  ffg,  respectively.  Their  concirrsion  regarding  these  three  desigrts,  based  oo 
the  trace  criterion,  rum  be  summarized  as  follows  (Muse  and  Anderson,  1978,  p,  162): 

(a)  !f  <7^  is  the  dominant  variance  component,  design  Dts  is  preferredi 

(b)  If  ai  ^ Og  with  either  ai  or  a}  larger  than  design  Dt^  is  preferred; 

(c)  If  both  tr,  anti  tr|  are  larger  than  of,  design  Hoa  is  generally  best. 

Their  conclusion  of  favoring  certain  unbalanced  designs  over  the  balanced  design  is 
Quite  interesting.  However,  we  have,  In  general,  little  prior  hnowiedge  of  variance 
components  themselves  which  we  would  like  to  estimate.  Moreover,  their  asymptotic 
results  do  not  nrxeasarily  hold  for  small  samples.  Therefore,  it  is  neceseary  to  inves- 
tigate the  overall  quality  of  estimation  for  the  desigrts  in  a more  comprehensive  way. 
This  can  be  accomplished  by  uring  QDGs. 

The  estimation  capabQity  of  the  designs,  Dts.  and  using  ANOVA 
estimators  of  of  and  o|  can  now  be  compared  using  the  corresponding  QDGs  for 
IV«a  and  Wgg  in  (3.42)  and  (3.43),  resperttlvely. 


The  balanced  design,  Ora 

The  ANOVA  table  for  model  (3.32)  with  Dts  is 


Source  d.f.  Sum  of  Squares  Expected  Mean  Squares 
A 5 SSa  «»2  + <»J 


B 5 SSe  + 

Error  25 SSj ef 

’Ibcal(adj.]  3S  55r 


The  scaled  ANOVA  estimalora  of  and  Ug  are,  by  Corollary  3.3.1,  distributed  as 


The  balanced  disjoint  rectangles  design.  Dus 

The  ANOVA  table  for  model  (3.32)  uang  Dta-  when  Typo  I suras  of  squares  are 
considered,  is 


Type  I SS 

Expected  Mean  Squares 

A 17 

SSt  = vQfV 

2i7j  + lfe|  + of 

B 9 

S3B  = y‘Qay 

2aJ  + <rf 

Error  9 

SSs  = y'Ry 

<r? 

Total(adJ.)  33 

SSr 

Model  (3.33)  can  be  written  as 

SI  = ylw  + (/,s  ® lj)a  + (i,  ® I, » h)a  + e. 
The  varlance<ct}varlaiicc  matrix  of  si,  E in  (3.34},  is  of  the  form 
E = {/i,  a J-,)al  + (/,  ® Jj  ® /j)oJ  + Ix”], 
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and,  hann,  in  (5.46)  aod  £3  In  (3.47)  are 
respectively.  Also,  Qg  in  (3.44)  and  in  (3.45)  are 

e«  = 

respectively.  Ttierernre,  using  the  matrices  Q„,  Qg.  E„,  and  Es,  the  distributions  of 
Wgg  and  WgA  ibr  the  design  are  obtained  from  Tiicorrm  3.3.1. 

The  off’diagonai  design. 

Similarly,  for  D^,  Che  ANOVA  table  is 


Source  d.l 

Type  I SS 

Expected  hfean  Squares 

A n 

SSt  = vQaV 

2i7?  + j|ri  + <r? 

B 12 

S3b  = vQbV 

M + <^ 

Error  6 

SSb  = y'Ry 

cr> 

Totsl(adj.)  35 

. SSt 

The  model  is 


' = ttl«  + (Tis  ® l2)cr  + {ts'SK)0  + s, 


0 0 1 
10  0 
0 0 1 
1 0 0 

Z = (/,.  8 J,K  + (U^KK-y,  + 
&nd,  henre,  snd  Are 

Ee  = + + 

respectively.  Also, 


UsQg  the  matrices  Qg,  Qg,  £o,  and  Zg,  the  dlBCributions  of  W^a  sod  IV^a  for  tbe 
design  D^gg  am  also  obtained  from  Theorem  3.0.1. 

The  qimncUes  of  Wgg  and  WgA  for  designs  D^,  Dta,  and  D^gg  are  computed  in 
conjunction  with  Davies'  {1980)  algonthm.  With  each  design,  values  of  Q^{p)  and 
Q^{p)iit  (3.50)  and  (3.51),  respectively,  and  values  ofQ^(p)  andO^(p)  in  (3.52) 
and  (3.53),  respectively,  ace  obtained  for  each  of  several  values  of  p (0  < p < 1).  Tbe 


quantile  values  of  Woa  fot  selected  values  of  p 


Table  3.6.  Minimum  and  maximum 
for  designs  D»,  Da^^,  and  find 


optimisation  orpg(p,piitb)  and  tvith  respect  to  (piiPs)  were  carried  out 

Sfn.e.)  = {(bi.'h)l  0.5  <«,  <10.0,  0.5  < fs  < lO.O}.  (3.56) 

More  speciOcaJly,  for  a given  design  and  a given  p,  9n(p,Oi>09)  ^ computed  using 
the  following  values  of  (171, pi);  tj,  is  ifj  s 0.5(0.1]1  and  1(1)10.0.  The  maximum 
and  minimum  of  the  resulting  values  of  9p(p.r?i.rig)  for  WaA  and  9n(p.qi,qa)  for 
WgA'  which  provide  approximate  valtics  of  <?o!?(p).  OK!0>).  and  Q3y(p),  0o5(P). 
respectively,  are  given  in  Tables  3.6  and  3.7.  The  QDGe  for  the  designs  Ch,  Dm, 
and  Dm  are  constructed  on  the  basis  of  tiiese  values. 


(a)  DJi6  <8  DJxl2  |b)  O.bS  vs  D_od3  (cl  DJxC  vs  0_gd3 


Figure  3.7.  QDGs  of  H'oa  tor  the  three  designs,  £>u,  Ac.  end  D, 


la)D_bSvsDJxl2  (blCLMvsO.olS 


Fi^rb  3.8.  QDGs  of  WgA  for  rbe  three  deigns. 


A display  of  the  QDGs  For  and  For  rnoh  pair  oF  designs  are  pven  in 
Figures  3.7  and  3.3,  respectively.  The  following  observations  can  be  made  From  these 

(a)  The  variability  in  the  quantile  values  of  \V^  = ^ with  respect  to  (171,03)  For 

design  is  much  smalJer  Chan  For  designs  Dtg2  and  Don  |scc  Figure  3.7(a,b)j. 
The  variability  For  design  Dm  is  relatively  smaller  than  For  design  Dm  [see 
Figure  3.7(c)].  Thus,  we  can  get  very  stable  Type  I ANOVA  estimate  of  with 
Dm  than  with  Dm  or  Dm-  Also,  we  can  have  more  stable  Type  f ANOVA 
estimate  of  with  Dm  than  with  Dm- 

(b)  The  variability  in  the  quantile  values  of  Wg^  ~ ^ with  respect  to  [171, rp) 
is  larger  for  designs  Dm  and  Dm  than  for  design  Du-  Moreover,  negative 
estimates  are  more  Frequent  with  the  unbalanced  designs  than  with  Du  [see 
Figure  3.3(a,b)].  Tlius  a niore  stable  estimate  of  o|  is  obtained  with  Du  than 
with  Dm  or  Dm-  The  QDGs  For  designs  Dm  and  Dm  are  quite  similar, 
although  the  variability  for  design  Dm  is  slightly  smaller  than  that  For  design 

[see  Figure  3.3(c)). 

(c)  The  degree  of  uniformi^  in  the  quantile  values  of  Wgg  is  larger  with  Du  than 

with  Dm  or  Dm  [see  Figure  3.3(a,b)|.  However,  there  is  no  noticeable  difference 
in  its  degree  of  uoiformlty  between  Dm  and  Dm  [see  Figure  3.8(c)). 

for  estimating  cither  or  Og-  The  unbalanced  design  Dm  is  slightly  preferred 
over  design  Dm  for  estimating  whereas  Dm  is  slightly  preferred  over  Dm  for 
estimating  o|. 


1.3  4 CoiinMriwii  of  dwitiM  uaiin  ML 


In  tiiie  section,  we  consider  the  probtem  of  design  comp&rison  for  model  (3.33) 
neing  ML  estimation  of  or  tr|.  Let  OgA/  ^6\t  denote  the  ML  estimotes  of 
or  o|,  respectively.  Let  us  also  denote  the  scaled  ML  estimates  of  and 
by  Wou  — and  Wgu  = respectively.  Note  that  and  do  not 
have  a closed'fonn  expression  regardless  of  whether  the  data  set  is  balanced  or  not. 
Furthermore,  the  exact  distribution  of  ML  estimatoia  is  unknown.  Theiefbre,  as  in 
Section  3.2.4,  we  need  to  use  the  EQDGs  for  design  comparisons.  An  outline  of  the 
corresponding  Monte<Carlo  simulation  folbws: 

(a)  Without  any  loss  of  generally,  we  assume  that  = 0 in  model  (3.32).  We  then 

Stot  = “.  + ^'r  + ‘u*.  i = l,2,— ,r,  j = l,2,.-,s,  i = 0,:,2,-,r^,. 

(3.S7) 

Fbr  the  balanced  dmign,  Dtg,  r = s = 6 and  ^ = I.  Fbr  the  unbalanced  designs. 
Duj  and  D^,  s D or  1 satisfying  N = = 36.  We  recall 

that  'll  = |!  and  ih  = ^.  Hence,  a,-  ~ 1^(0, mof),  A ~ 1^(0,  tsoj),  and 
c>>k  1V(0,^).  Let  us  denote  a region  in  the  parameter  space  of  the  triple 
(bii'U.i^)  by  V. 

(b)  For  a given  (ih,rb>i^)  la  V and  a given  design  D,  a random  vector  y is  generated 

on  the  basis  of  model  (3.37):  one  random  number,  which  we  denote  by  a,,  is 
generated  from  ^(D,<|i£T^}.  Independently  from  Ot,  another  random  number, 
which  wc  denote  by  6j,  is  generated  from  lV(0,'7s^).  In  addition,  one  random 
number,  wiiith  we  denote  by  Spi,,  is  generated  from  IV (0,  a?)  independently  of  o,- 
and  bj.  This  process  produces  an  observed  Pu*  value,  namely  pgs  = o,+6^  + su* 
which  represents  a response  value  in  the  t*^  row  and  they^  column.  For  designs 


Dt4i  &nd  the  generatiorL  of  Siji,  ie  made  only  for  the  (i,  j)^  cell  of  the 
incidence  matrix  X whose  element  correaponde  to  1 (i.e.,  rt(j=l).  The  incidence 
matrices  for  Dt^  and  were  tlefined  in  Section  3.3.3.  By  repenting  the 
process  of  generating  S of  with  appropriate  Oi's  and  hj's,  we  obtain  a 

(c)  UsiDg  the  response  vector  y obtained  in  step  (b),  ML  estimates,  and 

are  computed.  Values  of  and  Wgu  — are  then  obtained. 

(d)  By  repeating  steps  (b)  and  (c)  a sudicient  number  of  times,  we  can  obtain  empir- 

icai  quantiles  of  and  for  selected  values  of  p for  a given  (yi , ps,  of )• 
We  denote  the  p^  empirical  quantile  thus  obtained  by  ^(p,7i,Th,of)  and 
9p(p,hiirbiOf)  for  Ifou  and  tV^v,  respectively, 

(e)  Steps  (b),  (c),  and  (d)  ore  repeated  several  times  corresponding  to  different 
values  of  (rti,q3,of)  selected  from  the  region  V.  The  set  of  such  selected  triple 
is  denoted  by  T".  Thus  c V, 

(f)  The  maximum  aod  minimum  of  the  p**  empirical  quantile  vajuee  over  the  set  V 

in  step  (b)  are  obtained  for  each  of  fVqsr  and  Wgu~  Wc  then  have,  for  WqAr, 

of),  (3.58) 

QSS(P)  = , min  ffl(p,ih,'h,of),  (3.59) 

* (iheie?irr* 


(si»e?)e7- 

0S5(p]  = 


(3.81) 
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(g)  Values  of  Q^{p)  and  Q^(p)  are  plotted  against  p to  obtain  the  EQDGs  tor 
Wqu-  Also  the  values  of  Qoj  (>>}  and  Q^(p)  arc  plotted  against  p to  obtain 
the  EQDGs  for 

Let  us  novr  illustrate  the  application  of  the  above  procedure  using  the  desgns, 
Dies,  and  Dgus.  The  region  V is  chosen  as 

''“{(li.'fe.of)  lOSSbi  <10,  0.5<i)a<  10,  0.1  < of  <10} - 

(3.62) 

The  set  T*  mentioned  in  step  (e)  consisu  of  1125  points  selected  from  V sutb  that 
T‘ = {(th.ih.o,*)  I 1)1  = 0.6(0.1)1  and  1(1)10,  i)s  = 0.5(0.1)1  and  1(1)10, 
o?  = 0.1,0.5,l,5,lB). 

ASAS  Macro  (SAS,  1090)  was  written  to  compute  ^(p,f)j,  and^(p,t)i,t)s,^} 

values.  The  PROC  MIXED  procedure  in  SAS  (1997)  wee  used  to  obtain  ML  estimates 
of  Oq  and  in  the  SAS  Macro  program.  A total  of  2000  replications  in  step  (c)  were 
made  to  calculate  emplcicat  quantile  values  for  each  design  and  each  triple  (91,93,0^) 
in  the  region  TV  Table  3.8  gives  values  of  QS?(p)  »"'i  for  »»«•  Table  3.9 

gives  values  of  Qoj(p)  and  0o5(p)  for  Wgu-  EQDGs  plots  of  are  given  in 
Figure  3.9  and  those  of  Wgu  are  given  in  Figure  3.10.  We  note  that  the  variability 
in  the  quantile  values  for  the  balanced  design  Dm  is  smaller  than  those  for  designs 
Dm3  and  [sec  Figures  3.9(a,b]  and  3.10(a,b)].  However,  Che  distinction  between 
design  Du  and  design  Du  is  not  so  much  noticeable,  even  though  design  Du  is 
slightly  better  than  design  Du  >n  its  variability  [see  Figures  3.9(c)  and  3.10(c)).  These 
observations  are  consistent  with  those  made  with  ANOVA  estimates  in  Section  3.3.3 
Nate  that  Figures  3.9  and  3.19  appear  to  be  identical.  It  is  because  r = s for  all 
designs;  r s s = 6 for  Die,  end  r = e = 18  for  Du  end  Du.  Also,  in  the  region 
V in  (3.62),  ranges  of  91  and  93  are  identical,  so  that  the  EQDGs  for  Wau  must  be 
similar  to  those  for  U ssr  for  all  designs. 


Figure  3^.  eQDGs  of  I 


for  the  three  deigns,  Du,  D^,  uid  De, 
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(>)  OJlS  vs  OJld2  (b)  OJ>6  vs  D_Od3  (c)  DJxC  vs  D_I>I3 


Figure  3.16.  EQDGs  of  Wgn  for  the  three  deigns,  D»,  Dtgi,  and  Daji 

3.3.5  Comnafuoit  of  ANOVA  tod  Ml  ealmalioa 

We  can  compare  ANOVA  estimation  with  that  of  ML  with  regard  to  or  ir|. 
This  comparison,  however,  is  limited  since  the  regions  5(n,.n,]  In  (3.56)  and  V in  (3.62) 
are  not  the  same.  The  QDGs  for  and  the  EQDGs  for  Wou  are  shown  together 
in  Figure  3-11  for  each  design.  The  degrees  of  uniformity  In  the  quantile  values  of 

for  the  unbalanced  demgns,  Ossn  and  the  variability  in  the  quantile  values  of 
H'o*,  is  much  smaller  lhan  that  of  (see  Figure  3-U(b,c)].  The  QDGs  for  WgA 


find  thfi  EQDGfi  for  fino  6bo«ii  in  Figure  3.12.  It  is  interenting  Chnt  the  degrees 
of  uniformity  in  the  quantile  values  for  the  scaled  ANOVA  esciniates  find  for  the 
scaled  ML  estimates  are  very  Bimilar  for  all  deagns  over  the  range  of  p.  Moreover, 
the  variability  In  the  quantile  values  with  respect  to  (qi)  th)  ^er  the  scaled  ANOVA 
estiinfites  is  smaller  than  that  for  the  scaled  ML  estimates,  except  for  the  possibility  of 

In  summary,  the  ML  estimation  provides  much  more  stable  and  precise  estimate  of 

the  ANOVA  estimation  produces  more  precise  estimate  than  that  of  ML  regardless 


(a)  D_ba  (b)  DJXI2  lei  D.ecQ 


Figure  3.11.  QDGs  of  IVqx  and  EQDGs  of  W„u  for  the  three  designs,  Das,  D«c.  and 
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I EQDGs  of  the  three  de&^s,  Dte>  Dt^,  and 


3.4  The  Thrw-Pold  Nesttd  RMidom  Model 


staggered  i 


dosisn  spreads  the  degrees  of  Creedom  aliaost  evealy  over  alJ  the  stages.  Tliey  also 
mentioned  that  this  design  is  not  only  simple  to  administer  but  is  also  "open-ended", 
which  means  "the  design  allowa  more  levels  of  the  top  stage  variable  to  be  added 
without  upsetting  the  analyris"  (Leone  et  al.,  1966,  p.  737).  It  is,  however,  still  of 
interest  to  know  how  good  the  staggered  nested  design  is  as  compared  to  the  balanced 
nested  design  or  the  inverted  nested  design. 

The  estimation  problem  of  the  variance  components  for  the  staggered  nested  de- 
sign has  been  addressed  by  Babbridge  (1963,  1966),  Smith  and  Beverly  (1961),  Nel- 
son (1983.  1993a,  1993b),  and  more  recently  by  Khattree.  Nallt,  and  Mason  (1997). 
Nelson  (1993a)  asserts  that  the  staggered  nested  de^gn  is  t)ic  most  popular  among 
all  nested  designs.  Khattree  and  Naik  (1993)  presented  statistical  tests  for  variance 
components  m the  staggered  nested  derign. 

Figure  3.13  illustrates  part  of  the  four-stage  balanced  nested  design,  the  four-stage 
staggered  nested  design,  and  the  four-stage  inverted  nested  design.  In  this  section, 
these  three  designs  will  be  compared  in  a mote  comprehensive  way,  namely  by  loukbg 
at  the  disperwon  of  the  quantiles  of  the  variance  components  estimators.  The  deigns 

designs  will  be  shown  later  b their  ANOVA  tables. 

Consider  the  unbalanced  three-fold  nested  random  mode], 

Pyu  = jr  + a,  i-  ^(,)  +■  7k„)  -h  f„-u.  (3.63) 

^■1  ^*1  ^ti  = N.  Let  us  denote  the  factors  of  first 

three  stages  aa  A,  B(A),  and  C(AB),  respectively.  In  model  (3.63),  p is  an  overall 
mean,  a,  is  the  effect  due  to  the  level  of  factor  A,  is  the  ellect  due  to  they^ 


(a) 


^ jb) 


Figure  3-13.  Four-Stage  Nested  Designs:  (a)  the  BaUneed  Nested  Design  Showing 
The  Piist  and  Last  L^ls  of  the  Stage  1,  (b)  the  Staggered  Nested  De^gn  Showing 
The  Flist  and  Last  Leveis  of  the  Stage  1,  and  (e]  the  Inverted  Nested  Design  Showing 
The  First  Four  Levels  of  the  Stage  1 lAdapled  from  Leone  et  ai.  (1968)} 


w1w-^X,q  + X,(3  + X-j7  + «. 


(3.64) 


X,  = 1.,,. teeny 

X,  = 9j..®?..®2i.l^., 


S = W.(.),A(.,.-.  A, 
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and  ML,  where  dj,  and  denote  an  estimator  of  oj,  and 
respectively.  Three  designs,  namely,  a foiir-stagit  iwlanred  nested  dengn,  a four-stage 
staggered  nested  design,  and  a four-stage  inverted  nested  desgn,  are  consideretl.  The 
total  number  of  observations  for  eatdi  design  ia  *10.  Tiiese  desigua  are  the  ones  used 
by  Leone  and  Nelson  (1966)  and  Leone  et  al.  (196S),  so  that  we  can  compare  our 
results  with  theirs.  ANOVA  estlmatois  are  considered  in  Sections  3.4.3  and  3.4.4, 
whereas  Section  3.4.3  is  focused  ou  ML  estunators.  Contpariaone  of  these  estimation 
methods  fur  each  design  are  made  m Section  3.4,6. 

3.4.3  Quanlileaof^,  and  ^ uaina  ANOVA  estimation 

An  ANOVA  table  for  model  (3.63)  is  of  the  form 

Source  d.f.  Sum  of  Squares  Expected  Mean  Stiuares 
A 0-1  SSa  + 

B(A)  6.-0  SSs(a) 

C(AB)  c.. -6  + 

Error  jV-c,  5Ss  of 

Toial(a.U.)  N-1 Sft; 

S5a  = VQaV' 

SSb(a)  = v'QaV’ 

SSctAai  = vQcV' 

SSe  — y Rp, 


*■  = [3.72) 


by  and  rtspcciively.  Using  Theo- 


rem 3.2.1  and  Cornllftry  3.2.1,  we  oin  obroin  the  distribution  of  li'ax,  iV^^aiH.  ’tnd 

ThiHjrem  3.4.1  For  the  unbalanced  three-fold  nested  random  model  (3.63),  Ihe  scaled 
ANOVA  rrstimators  of  a\,  , and  distributed  os 


/ = number  of  distinct  positive  ^geiivalires  of 

g = number  of  distinct  positive  eigenvalues  of  Q^Sa, 

A = number  of  distinct  positive  ^grnvairtes  of  Q^S.,, 

7,  = positive  eigenvalues  of  <J„E„  with  multiplicity  I,,  j = 1,2,.-- ,/, 
ij  — j**  positive  eigenvalues  of  QaEa  with  multiplicity  mj,  j w 1,2,. ..  ,g, 
positive  eigenvalues  of  (J,E,  with  multiplicity  ’V,  j = 1,2, . . . .A, 


(3.76) 


(3.77) 


(3.78) 


Q. 


Q, 

<3, 


E, 


Qa- 


- k,kt) 


ik,k,  + k.k,-k,k,-k,k,) 

*,*<k.(,V-c)  “■ 

I ^ __*5_o  (t.-*.)  B 


75 


(3.79) 


(3.81) 


(3.82) 

= XiX'i  j + XjXj  + X,Xi  + /»  , (3.83) 

* XiX;  + XjXi  j + XjXi  + /„  .(3.84) 


Korf,  Che  ntHCrices  Q,,.  Q^,  Qc.  anti  Race  dehned  in  (3.66)  through  (3.69),  and  k^'s. 
I = 1. 2,'  ••  ,6,  are  the  cueflicienu  of  the  variance  componeota  In  the  expected  mean 
stiuares,  wliicli  are  defbed  in  (3.70)  through  (3.75). 

Proof  From  the  expected  mean  aquaret  pven  in  the  ANOVA  table,  the  ANOVA 
exCimatnra  of  ta  of  = A/Sg  = p'  (jj^)  V*  The  ANOVA  eatlmator  of  o|f„). 


= ;^IM5e,xa,-AfSil 

= " L(c  -5.)‘^'^"i.(X-c)H‘' 

= v'0,v. 


■ given  by 


= k!^{N-l.A  >' 


_ jk,k>*k,l^-k,U-kM 


loiwver  Q„1.V  = = Q,l„  = 0.  By  letting  X 

:aled  ,«OVA  eettaatore  of  a],,,,  and  a. 


:,  X),  = ^x:,and 


Then,  IV,.  = and  IV„.„.  = ^ ate  dl 


respectively,  where  9i  = ^,  % = and  jjg  — 

Proof  With  model  (3,d&),  the  expected  mean  squares  of  the  factors  are 
E{MSt)  = bauj’  + cnoJi,]  * + ‘'f • E(MSbia>)  = cno}(Bi  + ’“?(•»)  + ’’f' 

£(WSc(xa))  = tio^(ofl)  + “tid  £r{A/5ff)  — of.  Hence,  the  sums  of  squares  of 
the  faolocs  arc  distributed  as  SS<  - £(A/Sj)y2.„  SSbia,  ~ £(A/Ss(x))xftt-,|. 
^Sc(A0i  ■“  •^(AfJcxAO))  **  £{Af5s)  Vy,(n  ly  Moreover,  they  are 

Independently  distributed.  Therefore,  the  ANOVA  estimators  of  variance  components 


= j^IA/Sx  - iWSflinil 

0 I [E(MSa)  2 ^(MSbia,)  2 1 

- ten  [ 0-1  *-  0(6-1)  **'H 


Sines  li'oA  = H «(,,)/<  “ "*“l*  follows.  O 

Tbs  exact  distributions  of  W^a,  I^s(«)X*  ond  1Vi,(aS)4  >n  Theorem  3.4.1  end  in 
Corollary  3.4.1  ere  deteTmined  by  their  qiienriles,  whose  exact  vnlues  are  obtained  by 
uring,  for  example,  Davies'  (1980)  algorithm. 


The  quantiles  of  lV«a.  Ifst*)s<  and  depend  on  the  unknown  values  of 

hi  = * "^1  and  ift  = as  well  as  on  the  design  employed.  We  shall 

denote  a design  by  D.  Let  ua  denote  r;  = (q,,i7j,  nj).  The  dependency  of  the  p“ 
quantile  on  these  unknown  values  and  the  design  is  indicated  by  writing  9o(p.rf)  for 
If  aXr  ho(p.ri)  for  Wstaui.  and  9p(p,ri)  for  fkl^asia.  Let  us  also  denote  a region  in  the 
pHfainerer  space  of  ij  by  Sj|.  The  QDGs  of  W,a,  and  Wjjaa),,  are  obtained 

by  plotting  their  maximum  and  minimum  quantiles  against  p,  which  are  defined  as 


For  H'„. 


I3ob(p)  * ^“hSCP’r?))  and  (3J9) 

CSSOt)  = (3-90) 


For  Vi'sio)^. 


e~(P)  = »nd  (3.91) 

CSj(p)  = mhi  lilp.v)-  (3.92) 


C^tf)  = 9j(p,i|),  and  (3.93) 

Qn'W  = 

Aa  it  ia  stated  in  Section  3.4.1,  three  designs,  namely  a fdur*5tage  balanced  nested 
deeign,  a four-stage  staggered  nested  de^,  and  a four-stage  inverted  nested  design, 
are  to  be  compared.  We  denote  these  designs  by  Dm,  D,„,  and  D,„„  respectively. 
The  estimation  capability  of  tliese  designs,  using  ANOVA  estimation,  with  regard  to 
^o}<  ^tid  compared  using  the  QDGs. 

The  balanced  nested  design,  Dm 

The  ANOVA  table  for  model  (3.85)  using  design  Dm  with  a s 5,  5 — 2,  c = 2,  and 


Source  d.f.  Sum  of  Squares  Expected  Mean  Squares 
v'QaV  St'S  + ‘•"sta)  + 

B(A)  5 v’QnV  4o5,„,  + 2<rJ|,(,  + irj' 

C(.AB)  10  vQcV  2oJ,^-t-of 

Error  20  V Rv  o? 

Toml(adj.)  39 
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From  Corolloiy  3.4.1,  the  scaled  ANOVA  estimalore  for  and  are 

disiribuied  as, 

H'.a  = ^[8  + 4(^)+2(|)4-i]  xJ  (3.95) 


The  staggered  nested  design,  Dtig 

The  ANOVA  table  for  model  (3.83}  imng  design  D,tf  with  a — 10. 5.  w 20.  c.  w 30, 
and  /f  = 40.  is 


Source  d.f. 

Sum  of  Squares 

Expected  Mean  Squares 

A 9 

vQ.iV 

4‘’J  + H.)  + K»S)  + '^ 

B(A)  10 

vQaV 

C(AB)  10 

v'Qcv 

K.s)+'^ 

Error  10 

s ' fly 

of 

IdtaKadj.)  39 

SSr 

Model  (3.64).  in  vector  form,  for  riengn  Dag  can  be  written  as 

It  = /il«  + (Fit  ® li)o  + (Fio  ® -4)13+  (Jio®  B)“r  + e, 

where  A=  la  ®/i  and  B = Ij  9 Ig.  The  variance-covariance  matrix  of  it,  S,  is 
S = (/lo®  J()Oo  + (/lo « AAVIio)  + (Fio  ® BB')tr?,^a,  +/mo,’. 


c,  and  in  (3.82)  cbrough  (3.84) 


Sa  = (/,o®J.)(|)  + (/io®AA')  + (/io8Bfl')(Sj+/«(^^,and 
E,  = (/,o®J.)^^)  + (/io«AA')^^)*(/io®BB’)  + /«(^), 

respectively.  The  matrices  of  Q^,  Qfl,  and  Q,  in  equations  (3.79)  through  (3.81),  for 
this  d»ign,  are 

Qi  = 


Qa  = jail®  Jr)-;^J«, 

Qfl  = /io8X'- j(/,o®J.), 

Qr  = /io®B' -/loSX', 

B = /«-/io®B'. 

and  A'  = IJ3  ^ /i  and  B'  = ^ Ji  8 13.  The  matrices,  Q^,  Qy,  Q^-,  and  B,  are 
Dhtained  liy  using  equations  (3.66)  through  (3.69)  with  design  D,,g.  By  applying 
Theorem  3.4.1  vrith  the  matrices,  Q„,  Qg,  Q,,  E„,  S3,  and  E,,  the  distribution  of 
ll'aA,  iisfaiA,  and  for  the  design  B,,,  arc  obtained. 


N==iO. 


re  C = 1.  ffl  Is  © I2  © li,  D = li  9 li  e Ij  © X),  ami  £ = Ij  e Js-  1 


S = {I,®CC)oi  + (/.  ® DD')»J„,  4-  (/s  ® 


E„  Es,  and  E,  in  (3 


E^  = (XsaCC')(^)+(/,®DD’l  + (/.»££’)(^)+/«(;j),and 

E,  = (/.8CC')(i)+(/.8DD'l(|)  + {/.®££')  + /„(l), 


S3 


On  = 
= 

o,  = 


17020* 


Qa  = I,»C‘  - —Ju,. 
Qh  = i.8£>'-/.®C, 
Qf  = 7,®£'-/4®0‘. 


atid  C = ffi  jJj  © jJaffi/j,  D'  = jJj  ©/j  © jJa  © and  B'  = jJj  © /». 
Again,  the  mathraa  Qn,  Qc,  and  R ar«  ohuuncd  by  using  equatuma  (3.66) 
through  (3.09)  uoth  tlaeign  Dir„.  By  applying  Theoruni  3.4.1  with  the  matrices,  Q^, 
and  distribution  ofU^aAr  l^s(o)4,  and  for  the  design 

The  quantiles  of  li'aa,  iV^ajA,  and  ti'^(oS)A  for  designs  Dw  ^iis<  and  are 
eoiiipured  ueing  Davies'  (1980)  algorithm.  For  each  dnign,  values  of  maxiinuin  and 
minimum  quantiles  of  the  acalcel  ANOVA  estimators,  which  are  dchned  In  (3.80) 
through  (3.94),  are  obtained  over  the  region  of  the  parameter  space,  namely  Sj),  for 
each  of  several  values  of  p.  The  region  ia  chosen  according  to  the  one  used  1^  Leone 
et  al.  (1968),  ao  that  we  can  compare  our  tindings  to  theirs.  It  is  given  by 

5Tj  = {(rh, 02.03)1  I<rh  <9,  1<02<9,  l<o><9).  (3.99) 

.More  spcdfically,  for  a givea  p,  a quantile  value  Op(p,  r;)  ia  computed  umog  the  follow- 
ing values  of  t7  = (oi,  oa,  02)1  Oi  — 7(1)9,  02  — 1(1)9,  oa  — 1(1)9.  The  maximum  and 


S6 


Table  3.12.  Minimum  and  iDaximuiii  quamile  values  of  for  selected  values  of 

p for  desisns  Du.  D„„  and  D.^, 


la)  0_bal  «9  D_81g  D_bsl  vs  D_jnv  (c)  D_s)g  vs  Djnv 


8 of  for  ilic  three  designe,  D^,  jD.^,  and  Ding 


Figure  3.14.  QDGs 


(s)  D_Sal  vs  D_slg 


(b)  0 J»l  v»  OJnv 


(e)D_5IQVlDJlIV 


Figure  3.1S.  QDGs  of  for  the  three  designs,  £>w,  D,„.  end  Di„ 


(a)  0_MI  vs  O.sig 


(c)  D_slp  vs  OJnv 


|b)  OJlal  VI  DJnv 


Figure  3.16.  QDGfi  of  for  three  dedgus.  Du~  . and  D„ 
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diffrrence  is  nol  apparenl.  Thin,  » sliglitly  more  swblo  ANOVA  estimaU!  of  oj 
can  be  obtained  with  D,„  than  with  Did  I**  Figure  3.14(b)]. 

(b)  The  variability  in  the  quantiie  vaiuea  of  with  respect  to  ij  for  demgn  Did 

is  smolier  than  that  for  designs  D,ti  and  Dme  [see  Figure  3.13(a.b)j.  Design 
D,i,  shows  simlinr  variaUlity  to  design  D^i  (see  Figure  S.15(c)l.  Hence,  a more 
precise  esEimaut  of  is  obtained  with  Did  than  with  either  Dpg  or  Ddi- 

(c)  The  variability  in  the  quantiie  values  of  with  respect  to  tj  is  larger 

with  designs  and  Dmi  than  with  design  Did  [see  FigEtre  3.16(a,b}].  The 
smaller  degree  of  variability  with  Duj  than  with  D,„  is  more  apparent  here  [see 
Figure  3.16(b)].  Also,  Ddg  provides  smaller  variabiUty  than  D^i  [see  Fig- 
ure 3.16(c)].  The  degree  of  uniformity  in  the  quantile  values  of  W-,ipB]A  is  ab*o 
larger  for  Dm  than  for  An.  [««  Figure  3.16(b)).  Hence,  a more  stable  estimate 
offlS  is  obtained  with  Dm  than  with  either  Ats  or  A,n..  Also,  a more  stable 
estimate  is  obtained  with  Daa  than  with  Dmw 

In  suramary.  the  use  of  designs  Ai>  or  D„,  is  slightly  preferred  over  demgn  D^i 
only  when  the  ANOVA  eatimaie  of  o|[  is  of  interest.  This  is  expected  «nce  Ddf  and 
D„,,  have  larger  degrees  of  freedom  for  factor  A than  Dm-  However,  the  amount  of 
gain  uidng  D„,  or  Di„  over  Dm  is  net  large  enough  to  Justify  losing  balancedness. 
If  it  is  more  importaot  for  an  experimenter  to  obtain  a stable  ANOVA  estimate  of 
of  as  compared  to  o|jj,j  and  o^jg^p  then  it  is  recommended  to  use  D,n„  rather  than 
Dm-  Otherwise,  the  balanced  design  is  preferred  for  estimating  ofj„,  or  The 
preference  of  D„,  over  Da,,  occurs  only  when  is  to  be  estimated.  Otherwise, 
there  is  no  distcict  difleicnee  between  the  two  designs. 


respectively.  Let  up  also  denote  the  corTesponding  scaled  ML  estimates 
by  and  • Tha  EQDGs  of  H'„w, 

H'etniAi.  and  for  each  de^  are  obtained  uKOg  the  following  method  of  the 

Monte<Carlo  simulation: 

(a)  Without  any  losta  of  generality,  we  aaaume  that  p = 0 in  model  (3.63).  We  then 

Vut  = “i  + 9jl0  + tHd)  + ‘iju. 

1 = l,2,---,o;  ; = i = 1,2,---,Ch,  and  I = We 

recall  that  ifi  = ^,  l|j  = and  ife  = Hence,  a,  - A'(0,t||oj), 

8j(,i  ~ W(0,7/i^),  y»(„)  ~ jV(0,tht7f),  and  e,,u  ~ lV(0,irf).  Let  R denote  a 
region  in  the  paramewr  apace  of  the  ciiiadniplc  (fj.oJ),  where  tj  = (th,ib,  >h). 

(b)  For  a given  in  if  and  a given  design  D,  a rtindom  vector  y la  generated  on 

the  basis  of  model  (3.100):  one  random  number,  which  we  denote  by  n,  is  gener- 
ated from  lV(0,i7i^).  Independently  from  ri.  b,  independent  random  numbers, 
which  we  denote  by  sh.S(3.<  * * .Sa,,  are  generated  from  W(0,the<).  Indepen- 
dently from  Cj  and  Sy,  cy  independent  random  numbers,  which  we  denote  by 
.<u»o,  are  generated  from  fVfOitjjoJ).  In  addition,  iiyi  independent 
random  oumbent,  which  we  denote  by  ^ro  generated 

from  iV(0,o’)  mdependently  from  r„  a,,,  and  t„s.  This  process  produces  a data 
vector,  (r,-i-sn-ffm-i-ti,in.ri*t*a,i+liii+tt,iia,  • • - ,r,+s*-flt*,ea,+0(t,e,^iv,6,-,^)  • 
which  represent  response  values  From  an  i“  first-stage  nesting  group  (i  = 
1,2,  •••  ,o).  By  repenting  this  process  independently  a times  and  combining 


(3.104) 


ojid,  for 

(•■“I 

(g)  The  meximufn  and  minimum  quantile  values  obtained  from  (3.1CI1)  to  (3.106) 
are  plotted  ogoinst  p to  obtain  the  EQDGe  of  Wan,  and  W',(afl)M. 

respectively. 

Lei  us  now  illustrate  the  application  of  the  above  procedure  using  designs.  Dm- 
and  D,bc.  The  re^on  R is  chosen  as 

fl  = j(n.rf)  I 1 <i?i<9,  l<ijj<9,  1 <ih<9,  0.1  <i7?<10}- 

(3.107) 

The  set  H meniioned  in  step  (e)  con^ta  of  3645  points  selected  from  R such  that 

« = {(*».<^)l’!i  = 1(1)9.  th  = 1(1)9,  >6  = 1(1)9.  o?  = 0.1, 0.5, 1,5, 10}- 

A SAS  Macro  (SAS,  1990)  was  written  to  compute  5S(p,  i).  tr?),  jS(l'.'».^).  “<1 
9o(p.r|.of)  values.  The  PROC  MIXED  procedure  in  SAS  (1997)  was  used  to  obtain 
ML  estimates  of  variance  components  in  the  SAS  Macro  program.  A total  of  2000 
replications  in  step  (c)  were  made  to  calculate  empirical  quantile  values  for  each 
design  and  each  quadruple  of  (ri.of)  in  the  region  /f.  Tables  3.13, 3.14,  and  3.15  give 
those  maxima  and  minima  which  are  defined  in  (3.101)  through  (3.106)  for  designs 
D,ig,  and  D,„e. 

L'aing  these  tables,  the  EQDGs  of  listatit.  ®ttd  are  plotted,  os 

shown  in  Eiguces  3.17,  3.1S,  and  3.19,  respectively.  In  each  figure,  comparisons  for 
each  pair  of  designs  are  made.  From  these  figures,  we  note  similar  results  as  for 
AN'OVA  estimation;  The  unbalanced  nested  designs,  £>«,  and  An«.  provide  better 
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ML  estimatn  than  the  batonced  nested  design,  Dja,  in  terms  of  smaller  variabil- 
ity in  the  quantile  values  for  IVom  [see  Figure  3. 17(a,b)].  However,  in  estiinatiug 
“Toiiaai-  the  superiority  of  Dtai  to  D,,,  and  C,„,  is  apparent  [see  Figure  3.18(a,b)  and 
Figure  3.19(a,b)].  Design  D„,  is  preferred  over  D,„  only  when  the  ML  estimation 
of  Is  ronsidered  [see  Figure  3.l9(c)|.  Otherwise,  there  is  no  clear  preference 
between  them  [see  Figures  3.17(c)  and  3.U(c)|. 

In  cunclu^n.  the  staggered  nested  design  does  not  provide  an  overall  gain  in 
terms  of  its  variability  in  the  quantile  values  for  the  scaled  estimates  in  comparison 
to  Che  balanced  nested  design. 


(a)  0 Jal  vs  O_sto  (bl  DJ>al  vs  DJnv  (el  D_sig  vs  DJnv 


Figure  3.17.  EQDGs  of  W„m  for  tlie  three  demgns,  Dw  and  As 


(a)  D_td  VI  D_«g 


(bj  Djm  vs  DJnv 


(e)  Ojlj  vs  OJnv 


Fipire  3.1S.  EQDGs  of  for  tlireo  deigns,  Cw,  D,„,  and  D„ 


sure  3.19.  EQDGs  of  for  the  three  designs,  Du.  D,„,  and  Dm 


3.4.6  CQim>ifijBnofANOVA*ndML 


Comparisons  oF  ANOVA  rsttmstion  with  that  of  ML  with  regard  too^, 

^ ipade  iti  FigTires  3.20,  3.21,  and  3.32.  Again,  nolo  that  this  comparison  is 
limltod  because  the  regions  Sri  in  {9.29)  and  R in  (3.107)  are  not  Identical.  Aside 
from  the  posibilicy  of  having  negative  estimates  of  the  variance  components,  the 
ANOVA  estimates  are  comparable  to  the  ML  estimates  in  tlie  quantile  values  of 
and  [see  Figures  3.20  and  321|.  This  holds  true  for  all  three  demgna  It  is  noted 
that  the  frequency  of  getting  negative  ANOVA  estimates  of  and  is  very  liigh. 


that  for  lVl,(aS)sr  with  respect  to  (t;,  r^),  except  for  the  possibility  of  gettiag  negative 
entire  range  of  p. 


The  QDGs  associated  with  ANOVA  estimation,  and  the  EQDOs  for  ML  estima' 
tion,  provide  a powerful  grophical  tool  for  comparing  designs  for  estimating  variance 
components.  The  dependency  of  the  deigns  on  the  unknown  values  of  the  variance 
components  was  circumventetl  by  the  consideration  of  dispersion  of  the  quantiles  of  a 
given  estimator  over  a certain  region  of  the  parameter  space.  This  makes  it  posidble 
to  compare  designs  uniformly  with  respect  to  the  unknown  parameteis. 

Wc  compare  balanced  deigns  wltli  several  unbalanced  designs  for  various  models, 
namely  the  one-way  random  model,  the  tw>way  cross  clas^Bcatlon  without  biter- 


Chapter  2,  several  studies  have  shown  that,  usmg 
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and  EQDGb  cf  1 


for  tbe  ihr 
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(el  D_UMUlencei9) 


(c)  DJnv(lm«ned) 


Figure  3.21.  QD(^  of  Wji.j., , 


of  Wg^aj/j  for  the  three  designs, 


Hi)  D_slg  (slagssfM) 


Figure  3.22.  QDGe  of  and  EQDGs  of  for  tbe  three  deagns,  Dui. 

Ob,,  and 


balanced  design  is  not  always  favorable  to  estimate  variance  components  over  UO' 
balanced  designs.  These  studies  were  based  cm  tim  iarge  sample  properties  of  ML 
estimator  (Muse  and  Anderson,  1978;  Muse,  Andemon.  and  Thitakamol,  1982}  or 
on  the  empirica!  distribution  of  ANOVA  estimator  (Leone  and  Nelson,  1966;  Leone 
et  al.,  1968).  Based  on  our  study  to  asesa  the  overall  quality  of  the  eatlmation  abil> 
ity  rhrougli  the  QDGs  or  the  EQDGs,  we  can  determine  which  unbalanced  design 
provides  better  estimates  than  the  balant^d  design.  By  a better  design,  it  is  meant 
H design  which  provides  smaller  variability  iit  the  quantile  values  of  the  scaled  esti- 
mates, or  which  provides  a larger  degree  of  uniformity  of  such  values.  The  QDGs  or 
the  EQDGs  allow  us  to  compare  designs  lor  each  of  the  variance  components.  For 
example,  as  we  see  in  Section  8.4,  the  choice  of  a nested  demgn  can  be  decided  V 
which  variance  component  is  of  most  interest  to  an  experimenter. 

The  comparison  of  ANOVA  estimation  with  ML  estimation  is  also  made.  Even 
though  this  comparison  is  lunited  since  the  regions  of  the  parameter  space  are  not 
Identical,  it  has  shown  that  the  ANOVA  estimation  provides  estimates  comparable 
to  ML  estimates  except  for  tbe  poaibility  of  having  negative  ANOVA  estimates. 
This  is  an  interesting  result  itecause  it  is  based  on  small  sample  properties  of  ML 
estimation,  rather  Chan  on  an  asymptotic  beliavior.  Only  for  some  unbalanced  designs 
for  particular  models,  for  example,  a design  for  the  ooe-way  random  model,  and 
riesigns  Due  and  for  the  two-way  crou  clasmBcatlon  random  model,  docs  ML 
provitie  sixably  bettor  estimates  than  ANOVA  (see  Figures  3J(c)  and  3.11(b,c)]. 


CHAPTER  4 

INVESTIGATION  OF  THE  EFFECT  OF  IMBALANCE 
4.1  Tht  Probtbilitv  of » Neptivt  ANOVA  Enimale 
4.1.1  Intrwiuction 

h is  commonly  known  tiint  che  nnolyBls  of  variance  (ANOVA)  method  of  esti- 
mation of  variance  components  doee  not  preclude  the  pcenbUity  of  having  negative 
estimates.  This  ptoblem  has  been  discuased  several  authore,  particularly  in  the 
case  of  the  one-way  random  model.  See,  for  example,  Kelly  and  Mathew  (1993, 
1994),  Leone  et  ol.  [1966],  Prabhalmran  and  Jain  (1987),  Singh  (1969),  Tbn  and 
Wong  (1978),  Verdooren  (1982),  and  Wang  (1967),  to  name  just  a few. 

The  usefulness  of  ao  ANOVA  estimator,  d^,  of  a variance  component,  o’,  can  be 
measured  by  the  size  of  the  probability  P(d^  < 0).  Of  particular  Interest  in  this 
concern  is  the  ability  to  study  the  behavior  of  P{a^  < 0),  and  to  identic  the  factors 
that  in6uence  such  a behavior.  Some  of  these  factors  include  the  type  of  design  used 
and  the  true  values  of  the  variance  components  for  the  model  under  consideration. 
The  task  of  undertaking  this  sort  of  investigation  can  be  quite  daunting  when  the 
data  set  is  unbalanced.  Fbr  example,  it  ia  very  rllJBcult  to  determine  the  pattern  of 
imlMlance  that  adversely  affects  the  value  of  P(d^  < 0)  because  of  the  dependency  of 
this  probability  rui  the  modurs  variance  components.  The  difficulty  is  compounded 
by  the  non-uoiqueoess  of  ANOVA  estimates  due  to  the  availability  of  several  ANOVA 
tables  for  the  sanre  unbalanced  date.  Fiirthermore,  the  sums  of  squares  io  an  ANOVA 
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table  are,  hi  general,  neither  independent  nor  have  the  chi*8qiiarednee5  property  an  is 
the  ease  with  balanced  data  under  the  umml  ANOVA  asaumptiona. 

The  purpose  of  Section  4.1  is  ui  study  the  behavior  of  P{5^  < 0)  by  modeling  its 
values  in  terms  of  factors  which  are  considered  to  influence  such  a behavior.  Geo* 
eraliaed  linear  models  (GLMs)  techniques  ate  used  for  this  purpose.  The  proposed 
modeling  approach  provides  a better  insight  into  the  effects  of  data  iiDbaiance  on  the 
quality  of  Ab'OVA  estimation  of  variance  componenta.  fa  Section  4.1.2,  the  probabil* 
icy  of  getting  a negative  ANOVA  estimate  with  one-way  random  model  la  formulated. 
A method  for  generating  dt^gne  having  a speciSed  degree  of  imbalance  for  the  one- 
way model  is  mentioned  in  Section  4.1.11.  In  Section  4.1.4,  this  probability  is  modeled 
empirically  using  generalired  linear  models  techniques.  The  prediction  capability  of 

4 12  The  orobabilitv  of  a aesaiive  ANOVA  estimator  for  the  one-way  model 

Consider  the  uubalanced  one-way  raudom  model  (S.1)  and  its  vector  form  in  (3.2). 
The  vnriance-covarUnce  rnatria  of  the  vector  of  observed  values,  which  is  denoted  hy 
£,  ie  ahown  in  (3.3].  With  model  (3.1),  the  probability  of  a negative  ANOVA  estimate 
of  o’,  which  we  denote  by  d*.  can  be  found  using  Corollary  3.2.3,  namely 

f (.:  <•)-»■  (rh  gv  - (1^  'i-  < ")  ■ I‘‘1 

wliere  A^,'  • • .A)  are  the  distinct  pceitive  eigenvalues  of  QS'  with  multiplicities 
• • ,mr,  respectively i the  matrix  Q was  defined  in  (3.7),  and  Che  matrix  £' 
is  defined  as 

S*  = Jn,  + -/w. 


(4.2) 


whrrc  p = of  the  incquollty  inside  the  probability  (4.1) 

Is  a linear  combinatian  of  independent  central  ebi^squared  variates,  llie  probability  of 
a negative  ANOVA  estimate  of  oj  can  be  computed  for  specified  values  of  p. 

Let  us  denote  a design  for  model  (3.1)  by  D — (ni,ne,<  • • ,ns).  Note  that  the 
values  of  A^'e  and  the  corresponding  mt’s,  i — 1,2,. ••  ,r,  arc  delermlned  by  Che 
design,  D,  and  p.  Moreover,  os  will  be  seen  in  the  next  section,  a design  D for  model 
(3.1)  can  be  generated  by  specifying  the  total  number  of  observations,  /V,  and  Che 
degree  of  imbalance,  d,  when  the  number  of  treatments,  k.  is  fixed.  The  degree  of 
iiikbaiance,  is  a quantity  that  measures  departure  of  a design  from  the  completely 
balanced  design,  which  will  be  defined  in  the  next  section. 

The  exact  value  of  (4.1)  can  be  obtained  by  using  a computer  algorithm  given 
by  Davies  (19S0)  as  was  mentioned  in  Section  3.2.2.  It  should  be  noted  that  Che  exact 
probabliily  value  bi  (4.1)  is  not  an  explicit  function  of  the  parameters  S,  and  p.  It 
can.  however,  be  modeled  empirically  in  terms  of  Af,  and  p.  Using  this  probability 
as  a response  variable  and  the  parameters  N,  and  p as  control  variables,  we  can 
establish  such  a relationship  in  order  to  study  the  effect  of  imbalance  on  P{al  < D). 
The  proposed  modeling  of  P(a^  < 0)  depends  on  a method  for  generating  designs 
having  a specified  degree  of  imhaJanre.  This  will  be  discussed  in  Che  next  section. 

4J  J Generstl.m  nf  denienx  with  a Specified  decee  of  imbalance 


A design,  that  is,  a set  of  n,'vaiues  for  model  (3,1),  can  be  generated  for  given 
values  of  jV,  and  h.  lb  measure  the  degree  of  imbalance,  Ahrens  ond  Pincue  (1961) 
introduced  two  measures  of  Imbninncc  for  model  (3.1).  One  of  them  is  given  by 


where  g < 4 < 1.  The  smaller  the  d,  the  larger  the  degree  of  imbalance.  The  value 
1 is  attained  wiien  the  data  set  is  balanced.  Kburi  (1996)  introduced  a method  for 


Ikble  4.1.  Number  of  designs  genented  for  each  combination  of  N,  and  p 


generating  designs,  = {ni.nc,  • • • , n«},  having  specified  values  of  N and  d.  This 
method  will  be  used  to  model  < 0)  ns  a resimnse  function  of  N,  d.  and  p. 

To  demonstrate  the  development  of  the  aforementioned  relationabip,  let  ua.  for 
example,  consider  thel  ft  = 5 in  model  (3.1).  Thus  D « {Tt|,ns.ti3,n4,ns}.  Several 
combinations  of  N.  d.  and  fl  arc  selected  according  to  a 3^  factorial  arranpnient. 
Tlie  selected  values  are;  = 25, 50, 75;  d = 0-32, 0.60, 0,85;  p = 0.2, 0.4, 0.8.  For 
each  aasigiunent  of  the  triple  [Ai,  d,  p),  several  designs  are  gccerated  oo  the  basis 
of  Khurl's  (1996)  method.  The  application  of  this  method  was  carried  out  using  a 
computer  program  written  In  the  S<PLUS  language.  This  program  is  included  in 
the  appendix.  Note  that  the  generation  of  designs  requires  only  the  specification  of 
N and  d.  The  generated  designs  along  with  the  value  of  p determine  the  nonzero 
eigenvalues  of  the  matrix  Q£‘,  and  hence  the  probabijity  in  (4.1).  Table  4.1  lists  the 
number  of  designs  generared  for  each  of  the  27  combinations  of  N,  d,  and  p Unequal 
numbers  of  designs  for  each  combination  of  the  levela  of  IV  and  d are  generated  since 
the  availability  of  possible  designs  varies  according  to  the  values  of  N and  d-  The 
total  number  of  designs  thus  obtained  is  303. 


Empirioil  mftdriinf  of  Pig?  < Dt 


For  e»ch  sssignme nt  of  (N,  t,  p)  in  Table  4.1.  several  deagns  are  generated.  Some 
of  these  designs  arc  listed  In  Table  4.2.  In  this  table,  denotes  the  actual  value 
of  (h  for  a ^generated  design.  For  each  design,  the  nonzero  eigenvalues  of  Q£' 
are  obtained,  where  Q and  S'  are  defined  in  (3.7)  and  (4-2),  respectively.  The 
corresponding  value  ofP(d^  <0j  In  (4.1)  is  conipuled  using  Davies' < 1980]  algorithm. 
The  nonzero  eigenvalues  of  QE'  and  the  values  of  P(dJ  < 0)  corresponding  to  some 
of  the  generated  designs  are  given  in  the  Ibblc  4.3.  Since  several  designs  are  generated 
for  each  selection  of  IV  and  d,  the  values  of  P{dg  < 0)  corrcepondiiig  to  the  triple 
(fV.d.p)  represent  replicated  response  values  with  N,  d,  aud  p acting  as  control 
variables.  Let  us  therefore  denote  a value  of  P(d^  < 0]  by  p. 

The  data  in  Table  4.3  can  now  be  used  to  model  the  probability  of  a negative 
0^  against  N,  i,  and  p.  Given  the  nature  of  the  response  y,  being  a probability,  it 
would  be  appropriate  to  use  generalized  luiear  models  for  this  purpose.  The  logistic 

not  an  estimated  probabOity  obtained  as  the  ratio  of  the  nuiober  of  successes  to  the 
number  of  trials. 


p.  and  hence  on  ur.n,  for  each  triple  (iV.d.p)  in  Table  4.2.  the  mean  ITmi  and  sample 
variance  s]|„.  of  w„  for  the  i”  triple  can  be  computed.  Plotting  tvm  against  s^,, 


(4.4) 


which  maps  the  Interval  (0,1)  onto  [O.oo). 


against  Srtu.  for  the  various  values  of  t,  t = 1,2,  • • • , 27,  may  suggest  a possible 


Table  4.2.  Generated  desigua  for  tbe  one-way  model  with  h-S  [i,  denotes  the  actual 
value  of  d for  0 d-generatcd  dsdgn) 


Table  4.3.  The  nonzero  eipnralues  of  QZ‘  and  the  exact  value  of  P{al  < 0)  obtained 


On  the  basis  of  the  results  in  Thble  4.3,  the  rooffideiit  of  detcnnination,  for 
a simple  linear  regression  model  without  an  intercept  using  and  s^  ia  0.72  for 
m s 1.  0.70  for  m s 2.  and  0.S3  for  m s 3,  whereas  that  of  usbig  Wm  and  in,  is 
0.81  for  m = 1,  0.02  for  m = 2,  and  0.96  for  m = 3.  This  bidicates  a strong  linear 
relationship  between  the  mean  and  standard  deviation  of  un,  and  therefore  m = 3 ia 
seiected.  However,  since  v is  usually  not  large,  values  of  wa  are,  for  the  moet  part, 
close  to  zero.  We  therefore  multiply  wa  by  10’.  Now  let  ts  be  defined  as 


(iS) 


With  this  transformation,  we  eon  assume  that  the  coefflcicot  of  variacion  of  m is  con- 
stant. That  is,  Ov  is  proportion  to  (o^  oc  ftw),  where  and  Ov  are  the  mean  and 
standard  deviation  of  w,  respectively.  In  this  case,  with  a logarithmic  transformation 
as  a varianca-stabilismg  transformation,  we  can  use  ordinary  least  squares  method 
of  estimailon  in  order  to  obtain  estimatea  of  parameters,  but  the  intercept  te  then 
biased  (McCullagh  and  Nelder,  1989,  pp.  283-286).  On  the  other  hand,  under  the 
log-link  function  to  achieve  linearity  between  a response  variable  and  control  vari- 
ables and  n quadratic  variance  function  to  describe  the  relationship  bcCweeo  mean 
and  variance  of  w,  we  can  use  an  iteratively  weighted  nun-linear  least  squares  estims- 
tion  method  in  order  to  obtain  estimates  of  parameters.  For  the  quadratic  variance 
function,  it  is  meant  that  the  variance  of  nr  detrends  only  on  the  mean  of  ur  througii 
a quadratic  variance  function,  V'(p,o),  such  that  tr^  = 6rF(/i„],  where  dr  is  a constant 
dispersion  parameter  and  V'(p«)  = This  estimation  method  can  be  used  with 
assuming  that  ur  has  gamma  distribution  (see  McCullagh  and  Nelder,  1989,  p.  286). 
Therefore,  it  is  now  fearibic  to  choose  the  gamma  distribution  as  a probability  dis- 
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In  order  to  model  ag&inst  N,  d,  and  p,  a Ibik  function  la  needed.  Thir  Function, 
whlcii  we  donoto  by  g,  relates  to  tbe  soumlled  linear  predictor,  f,  which  ie  repre- 
sented by  a linear  model  in  IV.  d,  and  p.  Althou^  the  suggested  oanoiucal  link  for 
the  gamma  distribution  is  the  inverse  link,  we  choose  instead  a logarithmic  link 
log(pv),  where  log(pw)  is  the  natural  logarithm  oFPm-  The  reason  for  this  is  that  (or 
out  example,  the  inverse  link  produces  infeasible  results  as  will  be  seen  later.  Thus, 
the  basic  components  needed  for  our  generalized  linear  modeling  scheme  are 

(a)  The  response  variable  w in  (4.5)  is  asumed  to  have  a gamma  distribution  with 

(b)  The  linear  predictor  is  expressed  as  ^ where  x - (o,  d,  p)’  with  v =■ 

log(^),  ^ is  a vector  of  unknown  parameters,  and  f‘{x)  is  a vector  function  of 

This  way,  for  hf  = 25,50.75  (the  selected  values  oF  N in  Table  4.1],  y will  have 
the  values  3d219,  3.912,  4.317  with  s range  of  1.1. 

P.  = e‘.  (4.6) 

Let  p(x)  denote  the  predicted  value  of  p = < 0)  at  at  = (y,  d,  pY  inside  the 

d.  and  p in  Tbble  4.1,  that  is, 

C = {(y,d,p)  |!og(25)<y  <log(75),  0.32  < d < 0.B5,  0.20  < p < 0.80}.  (4-7) 

and  (4.6)  we  have 


iogli5(i)|  = /'(ce)fl. 


71og(10)  + 31og[ji^]  = /'(x)3, 

where  d Is  the  ntaximum  likellbuod  estimate  of  0.  Hence, 

li(x)  = U . (4.8) 

In  order  to  meseure  the  goodness  ofdt  of  model  (4.8),  we  use  the  scaled  deviance, 
or  the  scaled  Pearson's  chl-squared  statistic  (see  McCuUagh  and  Nelder,  1989,  pp.  33> 
34;  SAS,  1097,  pp,  263-284).  The  deviance  function  is  defined  to  be  twice  the  dif- 
ference between  the  maximum  attainable  log-likelihood  and  that  attained  under  the 
fitted  model.  The  former  is  achieved  sdth  a model  that  baa  as  many  parameters  as 

the  maximum  likelihood  estimate  of  the  dispersion  parameter,  which  for  the  gamma 
distribution  is  equal  to  the  variance  divided  by  the  square  of  the  mean  (see  SAS,  1997, 
p.  279).  The  scaled  Peaison's  chl-squarcd  statistic  is  obtained  in  a similar  manuer. 
I'nder  certain  regularity  conditions,  the  scaled  versiuns  of  these  two  statistics  have 
the  asymptotic  chi-squared  distribution  with  degrees  of  freedom  equal  to  the  number 
of  oirservations  (303  in  our  example)  minus  the  number  of  parameters  in  the  fitted 
model  (see  SAS,  1997,  p.  283).  Note  that  the  deviance  function  increases  as  terms 
ace  deleted  from  the  model.  Because  of  this  feature,  the  deviance  is  used  Co  compare 
ncxtixl  models  rather  than  as  an  absolute  measure  of  goodness  of  fit. 

Several  model  forms  are  conadered  for  the  linear  predictor  ( - /'(x)0.  These 
models  are  listed  in  Thbic  4.4.  The  corresponding  dovianccs,  scaled  devlanccs,  and 

Table  4.4.  The  CBNMOD  procedure  in  SAS  (1997)  is  used  Co  fit  these  models  oo  the 
baffis  of  the  logarithmic  Unk  Function,  ( = log(p,„). 

From  Table  4.4  we  note  chat  the  complete  secood-degree  model  has  the  smallest 
scaled  deviance  value  of  308.3168.  This  model  is  therefore  chosen  to  represent  the 
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TSblc  4.5.  Thi!  stBiuiard  fiTora  of  paramMfr  eatlmates  and  goodncss-of-fit  staiistic* 
for  model  (4.9) 


linear  predictor  in  (4.6).  In  this  caae.  f'(x)0  haa  the  form 

f(x)d  = U.8239  + 2.7088t'  + 0.9157«  + 2.816p 
-3.1081/#  - 3.1073i/<i  - 9.5592#p 
-0.43791/^  + 6.7049#’  - S.0771p’.  (4.9) 

Tile  standard  errors  of  the  parameter  estlmales  in  (4.9),  and  the  associated  P-valoea 
of  the  corresponding  teat  statiatica  arc  ^vcn  in  Tabie  4.5.  We  note  that  all  the  terms 
in  the  model  are  highly  aigniScant  with  the  eaccption  of  the  coeHiclent  of  #. 

Using  model  (4.8),  with  f‘(x)0  given  ns  in  (43),  we  can  compute  the  eatimaled 
prcbabiUtiea  of  a negative  estimate  of  at  the  27  factoriai  combinatioiis  of  (jV,  #,  fi) 
in  Ihble  4.2.  These  probabilities  ore  given  in  Ibbie  4.6.  Since  there  are  several  repli. 
cations  of  V W P(o^  < 0)  at  each  factorial  combination,  the  maximum  dilfeience,  aa 
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well  ae  the  average  difference,  between  the  vaiuee  of  y In  Table  4.8  and  the  eorie< 
eponding  estimated  probability  at  each  (S,i,p)  is  computed  to  cheek  the  adequacy 
of  fit.  The  restilta  are  also  shown  in  Table  4.6.  We  note  that  values  of  the  maxinium 
difference  range  from  -U.0283  to  0.0204,  and  the  values  of  the  average  difference  range 
from  <0.0084  to  0.0141.  These  values  along  with  chose  of  the  scaled  deviance  and 
scaled  chi-squared  statistics  in  Table  4.5  indicate  that  model  <4.8),  with  f{x]0  given 
by  (4.9).  provides  a good  fit  to  the  data  in  Table  4.3.  It  should  be  noted  chat  uang 
the  canonical  link  function  for  the  gamma  distribution,  namely  the  inverse  function, 
instead  of  the  logarithmic  bnk,  produced  infeasible  results.  This  la  caused  by  nega- 
tive values  attained  by  flx)0  at  some  of  the  aforementioned  factorial  combinations 
in  Table  4.2,  which  is  not  pusidble  snee  with  an  inverse  link,  f‘(x)0  would  be  the 
reciprocal  of  I0'(v{*)/[1  - »(*)]}“. 

The  modeling  of  the  probability  of  a negative  estimate  of  using  (4.8)  in  com- 
bination uilh  (4.9)  provides  a better  insight  into  the  effects  of  N,  and  p on  the 
quality  of  ANOVA  rstimation  of  Contour  plots  of  y[s)  in  (4.8)  can  be  used  for 
this  purpose.  These  plots  are  shown  in  Figures  4.1,  4.2,  and  4.3.  We  can  clearly 
see  that  the  probability  of  a negative  estimate  decreases  with  increasbig  values  of  n, 
d.  and  p [recall  that  v = log(W)].  In  particular,  for  fixed  values  of  v and  p,  as  in 
Figure  4.2,  the  probability  increases  as  tbe  data  set  becomes  more  unbalanced,  that 
is,  as  d decreases.  Note,  however,  that  a d value  as  low  as  0.60  can  still  result  in 
relatively  small  negative  probability  values. 

4.1.5  The  use  of  model  (4.81  for  desutn  comparison 

Seacle,  Caaella,  and  McCitlloch  (1992,  pp.  75-76]  considered  two  unbalanced  de- 
signs for  the  one-way  random  model  with  W = 25  and  k — 5.  The  purpose  was  to 
study  the  effect  of  imbalance  on  the  variance  of  the  ANOVA  estimator  of  a^.  The 
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■=  log(iV) 


Figure  4.2.  Contour  plole  of  the  . 


.imated  P(dp  < D]  for  ■■wh  level  of  d ]o=log(Af]] 
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Figure  4.3.  Contour  plole  (if  the  estimated  P{&1  < 0]  for  each  level  of  p [i/=log(N)| 


two  designs  are  D„  and  D,j,  which  were  defined  in  (3.20)  and  (3.21),  respectively. 
The  corresponding  values  of  0 ore  dai  = 0.51  for  and  — 0.2S  for  Thus 
D,2  is  more  unbalanced  than  D^i-  Along  with  these  designs  let  us  also  conader  the 
balanced  design,  Dt  In  (3,23),  for  which  dt  = 1.  For  given  values  of  ij  - and 
hence  of  p = we  can  calculate  the  exact  value  of  < 0)  for  each  of  the 
above  mentioned  three  designs  using  Davies'  (I960)  algorithm.  The  correspoudiog 
estimates  of  P{d\  < 0)  are  obtained  by  using  model  (4.8)  along  with  model  (4.9). 
Tlie  results  are  given  in  Thble  4.7.  Note  that  for  ds  — 1,  dij  — 0,28,  and  some  of  the 
p values  Id  this  table,  the  use  of  model  (4.8)  amounts  to  an  extrapolation  since  such 
valixea  fall  outside  the  tegiou  C b (4.7).  Nevertheless,  the  exact  values  of  P(dl  < 0) 
are  quite  close  to  their  corresponding  predicted  values,  as  can  be  seen  from  Tbble  4.7, 
Thus  the  prediction  capability  of  model  (4.8)  is  quite  good,  even  at  pmnts  ouUudc 
the  region  C.  FVom  Table  4,7  we  note  that  the  probability  of  a negative  estimate  is 
smaller  with  D.i  than  with  D^.  The  values  for  Di  are  smaller  than  those  for  D,\ 
and  Du3,  as  expected. 

4.1.6  Concludina  remarks 

The  use  of  generalized  linear  models,  io  combination  with  a method  to  generate 
designs  havbg  a specified  degree  of  imbalance,  can  be  quite  elTecCive  in  studying 
the  behavior  of  < 0)  for  the  one-way  random  model  Exact  values  of  such  a 
probabibey  can  be  obtained  using  Davies'  (1980)  algorithm.  These  values  arc  used 
to  fit  a model  in  which  the  response  is  P{dl  < 0)  and  the  control  variables  arc  i/,  d, 
and  p.  where  p w log(Af),  Af  Is  the  total  number  of  obser\iilions,  d is  the  measure 
of  imbalance,  and  p = If  such  a model  is  deemed  adequate,  Chen  it  can  be 
used  to  determine  the  effect  of  changing  N,  d,  and  p on  P(dJ  < 0),  and  hence  on  the 
quality  of  ANOVA  esUmation. 
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l^hle  4.7.  Prediction  capnbiiity  of  mode]  (4.8)  for  designs  Dt.  D^i, 
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It  is  also  possible  to  extend  the  proposed  ixiethodoiogy  to  hjgher<order  models. 
Khuri  (19S7)  introduced  a general  technique  to  measure  Imbalance  for  any  nested 
or  crossed-classilicatton  models.  The  extension,  however,  requires  generating  designs 
with  specified  degrees  of  imbalance  for  such  models.  The  method  introduced  by 
Khuri  (1996)  applies  only  to  the  one-way  model.  Extending  this  method  to  crossrd- 
classification  models  is  straightforward.  Howeter,  for  nested  models,  or  models  with 
crossed  and  nested  efieers,  the  extension  can  be  more  involved. 

Methods  of  estimation  of  variance  components  other  than  that  of  ANOVA  can 
also  be  addressed  with  this  proposed  methodology,  fbr  example,  we  can  consder 
minimum  norm  quadratic  uiibiaaed  estiraatioo  (MINQUE),  which  may  produce  neg- 
ative eslimates  (see  Searle.  Casella,  and  McCulloch,  1992,  Cb.  U;  Hao  and  Kleffe, 
1966,  Ch.  5).  In  this  case,  Davies'  algorithm  cannot  he  used.  Instead,  probabilities 
of  negative  estimates  of  variance  components  may  be  obtained  1^  simulation. 

4.2  The  Power  of  an  F-leg  With  Homogeneous  Error  Variance 
4.2.1  Introduction 

The  statistical  tests  of  variance  components  (or  a balanced  model  have  nice  prop- 
erties. More  specifically,  tbe  sums  of  squares  in  the  ANOVA  table  are  independently 
distributed  as  scaled  chi-squared  variates  imder  the  aasumptions  of  normality,  homo- 
geneity  of  variances,  and  independence  of  the  model's  random  effects-  This  enables  us 
to  have  exact  tests  fbr  variance  components.  We  can  also  have  optimum  teats  like  uiii- 
Eormly  most  powerful  (DMP).  uniformly  most  powerful  unWased  (UMPU),  uniformly 
most  powerful  invariant  (UMPI),  and  uniformly  most  powerful  invariant  unbiased 
(UMPIU).  These  tests  always  coincide  with  the  ANOVA-based  exact  tests  whenever 
such  tests  exist.  Wlien  no  exact  lest  exists  lor  testing  a variance  component,  we  can 


always  construct  an  approximate  F<test  using  Sattmhwaice's  procedure  (Satterth> 
waite.  1941,  1940).  Khuri’a  (1995a,  1995b}  work  coneeriung  cbis  approximation  is 
noted.  We  also  would  like  to  mention  Seifert's  tests,  which  often  yield  exact  unbiased 
tests  when  exact  optimum  tests  do  not  exist  (Seifert,  1981). 

In  tire  case  of  an  unbalanced  model,  exact  tests  of  variance  components  are  avail- 
able for  a limited  number  of  models.  These  exact  F-tests  are  based  on  an  idea 
originally  due  Co  Wald  (1947),  and  are  extended  by  Seely  and  El-Bassiuujii  (1983). 
However,  the  existence  of  a uniformly  optimum  test  is  rare  for  the  unbalanced  model, 
and  hence  locally  optimum  tests  have  been  developed.  For  example,  in  the  case  of  an 
unbalanced  one-way  random  model,  there  is  no  unifbnnly  optimum  test  for  testing 
the  variance  cojnponem  of  a treatment  factor,  but  a locally  best  Invariant  (LBI)  test 
is  available  (Das  and  Sinha,  1967],  On  the  other  hand,  El-Bassiouni  and  Seely  (1996) 
proposed  a modihed  harmonic  mean  test  and  showed  by  comparing  empirical  powers 
that  it  outperforms  the  locally  most  powerful  (LMP)  test  and  Wald's  F-tast  when  the 
design  is  extremely  unbalanced.  Therefore  a study  of  the  power  of  a test  is  still  of 
interest.  For  more  details  concerning  exact  and  optimum  tests,  we  refer  to  Ferguson 
(1967),  Lehmann  (1986),  Kariya  and  Sinha  (1989),  and  Khuri,  Mathew,  and  Sinha 
(1998). 

Note  that  the  power  of  the  ANOVA-based  F-test  depends  on  the  design  used  and 
the  values  of  variance  components.  Donner  and  Koval  (1989)  studied  the  e9«t  of 
imbalance  on  the  power  of  the  F-test  for  the  one-way  random  model.  They  compared 
its  empirical  power  to  the  corresponding  likelihood  ratio  test  (LOT),  using  designs 
having  various  degrees  of  imbalance.  Tliey  showed  that  the  LOT  is  more  powerful 
when  the  design  is  extremely  unbalanced.  In  Section  4.2,  we  shcpw  that  the  exact  value 
of  the  power  of  the  F-test  can  be  obtained  by  using  Davies’  (1980)  algorithm  since  the 
power  can  be  expressed  os  a probability  of  linear  combinations  of  independent  chi- 
squared  variates.  Using  this  exact  value,  we  model  the  power  in  terms  of  the  design 
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and  the  variaDce  componeats  to  get  a better  inaigbt  of  the  effect  of  these  factors  on 
the  power.  By  this,  we  can  investigate  the  effect  of  imbaiancc  on  the  powet  of  the 
F-test  for  a variance  component.  Moreover,  the  proposed  niodeiing  technique  can  be 
used  effcctiveiy  when  an  experimenter  wants  to  find  a design  which  ensures  a certain 
level  of  the  power  of  an  F-test.  This  can  provide  a guideline  to  an  experimenter  to 
deagn  statistical  experiments,  in  Section  4.3.2,  the  power  of  the  F-test  for  a variance 
component  with  one-way  random  model  is  formulated.  This  power  value  is  modeled 
using  generaliaed  linear  model  techniques  in  Section  4.2.3.  To  perform  this  analysis, 
we  employed  the  same  generated  designs  which  were  used  in  Section  4.1.  Concluding 
remarks  are  given  in  Section  4.3.4. 

4.3.3  Power  of  the  F-test  for  ai  for  the  one-way  model 

Consider  the  unbalanced  one-way  random  model  (3.1)  and  its  vector  form  m (3.2). 
The  variance-covariance  matrix  of  v<  the  vector  of  observed  values,  which  is  denoted 
by  £.  is  shown  in  (3.3),  The  sums  of  squares  for  the  treatment  factor  and  the 
error  term  are  ^Qy  and  ^Ry,  where  (J  and  H are  defined  in  (3.7)  and  (3.8). 
respectively.  These  two  sums  of  squares  are  independent  and  ^\^Ry  has  a central 
X^-distribution  with  h degree  of  imedom.  Also,  under  ^y*Qy  follows 
a central  distribution  with  k - 1 degree  of  freedom.  This  gives  us  a test  statistic 
for  testing  the  hypothec  as 


which,  under  Hq,  has  the  F-distribution  with  k - 1 and  N - k degrees  of  freedom. 
This  test  is  L'MPI  If  the  data  set  is  balanced. 

As  was  mentioned  bi  Section  3.2,  the  treatment  sum  of  squares,  y'Qy,  can  be 
expressed  as  a linear  combination  of  independent  central  chi-squared  variates  of  the 


(4.10) 


(4.11) 


v'Qv  = j^>^xl^  1 

where  Ai.ii, •••  ,Ar  are  the  distinct  poMtive  eigenvalues  of  QL  with  multiplicities 
ni|.mv,  • • • .Ttir,  respectively.  UNiig  the  relationship  in  (4.11),  the  power  of  the  F- 
test  for  the  hypothesis  can  be  expressed  as 

= pfE-'rxiiu-Aie.  (^)  Xw-sXJIfJlioj  , (4.12) 

where  A)  = j = 1.2, . • - ,r,  which  are  the  distinct  positive  eigenvalues  of  QE' 
with  multiplicities  rtij,  i = 1,2,--  - ,r,  respectively,  such  that  = * - 1.  The 

matrix  £•  = ^E  is  shown  in  (4.2).  Also  A;*,  = and  p = The 

exact  value  of  (4.12)  can  be  obtained,  again  by  using  a computer  algorithm  given 
by  Davies  (1D80).  This  exact  value  can  be  modeled  in  terms  of  IV,  d,  and  p,  in  order 
to  acquire  a deeper  insight  into  the  effect  of  these  variables  on  the  power. 

4.2.3  The  Effect  of  Imbalance  On  the  Power  of  the  F-lest 


A design,  that  is  a set  of  n,--values  for  tlic  model  (3.1),  can  be  generated  for  given 
values  of  N,  d,  and  k using  Khuri's  (199G)  ineihud,  as  was  stated  in  Section  4.1.3. 
With  the  talue  of  p = and  the  generated  design  D,  we  can  compute  an  exact 
value  of  the  power  of  the /■-lest  in  (4.12).  Since  this  power  is  not  an  explicit  function 
of  the  parameter  N.  d,  and  p,  for  a given  k,  it  will  be  modeled  as  a response  function 
of  these  parameters. 
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Let  us,  for  example,  consider  tr  = 3,  so  chat  D = {ni,  02,03,04.  ne).  Consider  a 3^ 
factorial  arrangement  having  the  selected  levels  of  = 25, 50, 75,  i = 0.32, 0.60, 0.85, 
and  p = 0.2,0.4,0.S.  The  303  designs  generated  in  Section  4.1.  which  are  listed  in 
Thble  4.2,  are  used  in  this  study.  The  power  in  (4.12)  is  calculated  using  Davies’ 
(1980)  algorithm  for  eacli  of  the  generated  design.  In  this  study,  a = 0.05  is  used. 
Part  of  the  generated  designs  and  the  corresponding  powers  are  listed  in  Table  4.8, 
Since  several  designs  are  generated  for  each  specification  of  N and  we  have 
replications  of  the  power  value  for  each  combination  of  N,  i,  and  p.  The  data  in  Th- 
ble  4.8  can  now  be  used  to  model  the  power,  p,  in  terms  of  these  three  variahles,  N.  d. 
and  p.  Conddering  the  nature  of  the  response  p,  being  a probalnltty,  it  would  be  ap- 
propriate to  confer  a generalized  linear  model  6ccicg,  as  it  is  stated  bisection  4.1.4. 
Since  0 < d < 1 and  0 < p < 1,  we  consider  the  logarithnuc  cransfonoatiou  of  S, 
u = log{N).  instead  of  N itself.  The  range  of  the  variable  v is  about  1.1,  which  is 
comparable  to  the  ranges  of  d and  p.  The  region  of  experimentation,  C.  la  described 
in  (4.7).  Let  us  now  consider  the  transformatiou. 


which  maps  the  interval  (0,1]  onto  [0,cx2).  Since  we  have  replicates  of  p for  each 
combination  of  triple  (u.d.p),  the  sample  mean  and  the  sample  variance  of 
uTft,  for  the  i*^  tripie,  i = 1,2,  • • • , 27,  can  be  computed  for  varioue  value  of  m.  A 
simple  linear  regression  model  without  an  intercept  uring  iOru  and  or  wn,  and 
is  Sued  to  project  a telaiionshlp  between  the  mean  and  standard  deviation  of 
Wm  or  between  the  mean  and  the  variance  of  Wn,  respectively.  The  coefficient  of 
determination,  using  ulmj  and  e^,  is  0.55  for  m = 1,  and  0,63  for  m = 2;  whereas 
that  of  using  uin,  and  Sn,  is  0.74  for  m = 1,  and  0.90  for  m s 2.  This  indicates  that 
lliere  is  a strong  linear  relationship  between  the  mean  and  standard  deviation  of  w„ 
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with  m = 2.  Therefore,  with  the  tronsfonoation. 


the  gomma  distribution  can  be  chosen  as  a probability  distribution  for  the  random 
variable  w.  Let  us  denote  the  mean  and  variance  of  in  by  ttv  and  respectively.  We 
choose  ( - /'(z)/3as  the  linear  predictor,  where /3  is  a vector  of  unknown  parameters, 
and  f(x)  is  a vector  function  of  x such  that  f'{x)B  is  a polynomial  inx  = (v,d.p)of 
a certain  degree.  Tb  relate  the  linear  predictor  ( to  the  mean  response  /a,,  we  consider 
the  logarithmic  link  function  for  reasons  smilai  to  those  stated  in  Section  4.1.4.  We 
thus  have  the  following  components  for  our  model  fitting; 

(a)  For  the  random  component,  the  respoose  variable  to  in  (4.13)  has  independent 
gamma  distribution  with  a mean  and  a variance  function  V(/t)  — 

(h)  For  the  systematic  component,  the  linear  predictor  ^ is  given  tb'  C = /^(*)i3, 

(c)  The  link  function  wlikh  relates  to^  is  the  logarithmic  function,?  = log(ji,). 

I^r  the  linear  predictor,  ?,  we  conrider  the  second  degree  model, 
f{z)0  = 3a  + 0\x  + z'Bx  . 

where  /S',  = (ft, A, ft)  and  B = ((ftj}),t,j  = 1,2,3.  The  GENMOD  procedure 
in  SAS  (1997)  is  used  to  fit  various  model  forms.  The  scaled  dovinnee  and  the  scaled 
Pearson's  clu-squared  stattstic  of  the  fitted  mode!  can  be  used  os  measures  of  the 
goodness  of  fit.  Based  on  these  values  and  the  significance  of  Wald's  test  statistic  for 
each  parameter,  the  model  selected  to  represent  the  linear  predictor,  ?,  is 

f(x)0  = 11.3122-  1.3815i/-2.27184-  2.4343p 


- 1.688Ue«  - 1.523Svp  - S.4999«p 
+ 5.4755(1^  - 1.6730p''. 


(4.14) 


Table  4.9.  The  standard  crroia  of  parameter  eatimatea  and  goodness-of-St  statistics 
for  model  (4.14) 


Let  p(oi)  denote  the  eatimated  value  of  y in  (4.12)  at  a given  e inside  the  region 
of  experimentation  C in  (4.7).  Using  model  (4.14),  it  can  be  calculated,  at  a given 


• (4.15) 

These  estimated  values  of  p for  given  specifications  of  i are  listed  in  Table  4.10. 
Since  there  arc  several  replications  of  p at  each  factorial  combination,  the  maximum 
dllTerence,  as  well  as  the  average  difference,  between  the  power  value  ofp  in  Table  4.8 
and  the  corresponding  rstimated  power,  {/{x),  at  each  triple  (ir,  d,p)  is  calculated. 
The  results  are  shown  in  Thble  4.1D.  We  note  that  values  of  the  maximum  difference 
range  froni  -0.0269  to  0.0396,  and  values  of  the  average  difference  range  from  -0.0137 
lo  0.0062.  These  values  ore  very  small  as  compared  to  the  siae  of  the  actual  power 
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values.  These  values  and  the  gODdneas  of  fit  lest  statistics  in  Table  4.9  indicate 
model  (4.L5)  with  f\is)0  given  by  (4.14)  explains  the  data  in  Table  4.8  very  well. 

To  see  the  elfect  of  the  explanatory  variables  w on  the  poiver  of  the  F-t«st,  cnnlour 
plots  of  y(«)  are  displayed  in  Figures  4.4,  4.5,  and  4.6.  The  power  increases  as  we 
increase  values  of  v,  d,  and  p.  However,  for  a fixed  value  of  one  variable,  we  ttao 
increase  the  power  by  increasing  a value  of  another  variable  and  decreasing  a value  of 
the  other  variable.  For  example,  we  see  in  Figure  4.4(bl  that  when  r/  s log(50),  the 
estimsteil  power  is  0.5  when  p = O.S  and  d = 0.65,  and  it  reaches  0.9  when  p - 0.7 
and  d = 0.45.  Thus,  when  ^ = 50  (i.e.,  p s log(50)),  we  can  incrcaae  the  power 
from  0.5  to  0.0  by  increasing  p value  from  0.2  to  0.7  but  decreasing  « value  from 
0.65  to  0.45.  Fbr  more  practical  aspects,  we  can  have  the  following  observations  ho' 
exaniiiiing  Figures  4.6(b)  and  4.6(c):  when  p = 0.4,  in  order  to  get  power  value  of 


0.9,  we  need  a design  with  JV  = 75  (he.,  p = 4.3)  observations  whose  meesure  of 
imbaJance  la  almost  as  hig^  os  d s 0.8  [see  Figure  4.6(b)].  However,  when  p = 0.8, 
a design  with  N - ZO  (he.,  p = 3.4)  and  d = 9.4  is  sufficient  to  get  the  power  value 
9.9  [see  Figure  4.6(c)].  For  another  aspect,  if  ws  llhe  to  ensure  the  power  of  an  F-test 
to  be  at  least  0.9,  a design  having  a degree  of  Imbalance  es  low  as  d — 9.4  is  good 
enough  when  we  have  p = 0.8  and  iV  s 30  (i.e.,  p s 3.4)  [see  Figure  4.6(c)).  These 
observations  can  provide  very  practical  guideline  lo  an  experimenter  who  is  in  a stage 
of  designing  a statistical  ex]>enrnent. 

Now.  the  question  is  how  good  the  model  is  to  predict  the  power  of  the  F-lest 
at  points  outade  the  experimental  re^n  C in  (4.7).  To  answer  this  question,  we 
again  consider  three  designs,  D^,.  and  D^s,  which  were  used  in  Section  4.1.5 
for  similar  purposes  to  investigate  the  effect  of  imbalance  on  the  probability  of  a 
negaiive  A.\'OVA  estimator.  For  each  design,  the  exact  power  of  the  F-test  for  given 
values  of  0 s tg  p s = -^),  is  caiculated  using  Davies'  (1960)  algorithm. 


Also,  the  corresponding  estimated  power  y{x)  based  on  the  model  (4.15)  with  f'{x)B 


Figure  4A.  Contour  plots  oF  the  estimated  power  of  the  F-test  For  an  anioiig*group 
variance  component  for  each  level  of  t/  ^ log(N) 


1J5 


Figure  4.5.  Contour  ploM  of  the  eetinmted  power  of  the  F-test  for  an  among-group 
variance  component  for  each  level  of  0 [r^log(yV)] 
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Figure  4.6.  Contour  plots  of  the  retimated  power  of  the  F-test  for  an  among-group 
varianoc  component  for  each  level  of  p [o=log(^)j 
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given  by  {4.14}  is  calculated.  These  values  along  with  their  diHereiices  are  listed  in 
Table  4.U.  Note  that  ^ = i for  the  design  D^,  s D.23  for  the  design  and 
some  p values  (/>  s 0.9091,  p = 0.0909)  fall  out^e  the  region  Cln  (4.7).  Nevertheless, 
the  predicted  powers  with  model  (4.15)  arc  veiy  close  to  their  cortespoiidiiig  exact 
values.  Thus  the  prediction  capability  of  model  (4.15)  along  with  model  (4.14)  Is 
good,  even  at  points  outmde  the  region  C.  As  we  know,  for  a balanced  design,  the 
F'test  is  UMPI  and  therefore  the  powers  for  should  be  larger  than  those  for  D^i 
and  D^.  We  can  sec  this  from  the  Table  4,11.  We  also  note  that  the  power  is  larger 
with  Dal  than  with  Das  for  all  of  p values. 

4.2.4  Ctmcluding  remarics 

With  an  unbalanced  one-way  random  model,  the  ANOVA-based  F-test  for  testing 
the  significance  of  the  random  effect  variance  component  ia  valid.  However,  this  test 
does  not  have  any  optimum  properties.  There  is  also  a locally  optimum  test  which 
is  different  from  the  F-test.  However,  the  use  of  the  F-test  is  still  popular  due  to 
its  simplicity.  In  this  respect,  the  power  of  the  F-test  is  always  of  Interest  for  an 
experimenter. 

The  power  of  the  F-test  for  a varianen  component  with  an  one-way  raadom  model 
depends  not  only  on  the  design  used,  but  also  on  its  variaoce  componenU.  The  exact 
value  of  the  power  can  be  obtained  using,  for  example,  Davies'  (1980)  algorithm.  In 
Section  4.2,  the  power  is  modeled  as  a function  of  the  total  number  of  observations, 
the  degree  of  imbalaner.  and  the  ratio  of  the  variance  componcols.  The  fitted  model 
can  be  used  to  see  the  effect  of  these  variables  on  tbe  power,  and  hence  on  the 
ANOVA-based  F-test.  Moreover,  the  fitted  model  in  (4.15)  along  with  model  (4.14) 
can  provide  guidelines  for  choosing  a design  capable  of  producing  desirable  power 
values.  This  modeling  technique  can  therefore  be  used  effectively  at  the  demgnbg 
stage  of  a statistical  experiment. 
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Nlictioii  capablliCy 
6,6,5};  A'=»25(i' 

f model  (4.15)  for  dsignE  D 

= 3.219),  d,  = 1.0,  p = ii: 

.(-S)  . 

9(=)'  it  i(»)-v 

2 0 6667 

0.93541  0.99363  -0.00822 

1 0.5000 

04^  0527W  ^^5 

O.I  0.0609 

0.15098  0.14753  0,00344 

1,11,11);  jV  = 26( 

= 3.219),  = 0.6102,  p = 

"(-S)  - 

i(*l  > iM-« 

02^7  o’^  o'm2M 

1,1,21);  JV  = 25(i. 

= 3.219),  rt«  = 0.2809,  p = 

.(-0)  . 

fi(*]  V li(*)-9 

0.46640  040080  0.06560 

0.26  0.2000 

0.16193  0.13752  0.02441 

0.1  0.0906 

0.10382  0.08158  0.02224 

X)  » the  prcdicud  vtiue  ol  Lbe  power,  p. 

The  extonsioR  of  the  proposed  methodology  to  higher*ordcr  models  is  of  interest. 
Investigation  of  the  effert  of  imbalance  on  the  power  of  tests  other  than  the  F-test 
is  also  of  interest.  For  example,  the  performance  of  the  modified  harmonic  mean 
test  suggested  by  El-Bassiouni  and  Seely  (11196)  can  be  examined  more  thron^out 
in  terms  of  the  power  using  the  methodology  shown  in  this  section.  In  this  case, 
however,  we  cannot  use  the  Davies'  nlgnrithin  to  calculate  the  power  of  the  test,  but 
it  can  be  obtained  by  computer  simulation. 

4.3  The  Power  of.an  F-test  With  Hetorocneous  Error  Variances 

4.3,1  Introduction 

In  many  experimental  situations,  the  error  variances  in  a pven  random  modei 
may  be  heterogeneous  due  to  a variety  of  reasons.  Inequalities  of  such  variances  can 
seriously  affect  an  ANOVA  test,  particuiariy  when  the  data  set  under  consideration 
is  unbalanced  (see,  for  example.  Miller,  1986,  Section  3.7).  This  problem  has  been 
investigated  by  many  authors  for  Ibted-effecls  models  (see,  for  example,  Brown  and 
Forsythe,  1974:  Bishop  and  Dudewicz,  1978;  Levy.  1978a,  1978b;  Krutchkolf,  1988; 
Wcerahandi,  1995).  Fbwer  articles,  however,  considered  situations  involving  random 
or  mixed  models.  Rao,  Kaplan,  and  Cochran  (1981)  discussed  poiut  estimation  of 
the  parametets  of  the  one-way  random  model  with  unequal  error  variances.  Jeyarat- 
nam  and  Othman  (1985)  provided  an  approximate  test  for  the  hypothesis  Ha-ol—O 
using  the  Satterthewaite's  approximation  (o  the  distribution  of  the  ainong-group  sum 


effects  of  iieteruscednsticity  and  data  imbaiance  on  the  probability  of  a negative  es- 


thc  usual  F-test  concerning  the  hypoiheris  Power  values  were  computed 


model.  Singh  and  .loshi  (1989)  investigated  the 


tiinnte  of  Oq.  Singh  (1991)  studied  the  same  ejects  with  regard  to  tbe  power  of 


u^ng  a doubly  infinite  eerice  of  incomplete  beta  functiona.  He  concluded  that  het- 
eic«cedaatieity  increases  the  power  of  the  test  in  dtuations  where  the  larger  veuiances 
ore  associated  with  the  smaller  group  aiees.  If  the  larger  tnriances  are  associated  with 
tiie  larger  group  sizes,  tlie  power  is  decreased.  A small  increase  in  the  power,  due  to 
heteroscedastlcity,  was  noted  with  balanced  data. 

in  Section  4.3,  we  model  the  power  of  the  usual  F*test  concerning  empirically 

modeling  scheme  provides  sdded  insight  into  the  effects  ofimbatanec  and  heterogene- 
ity of  error  variances  on  the  power  of  the  T-test.  In  particular,  the  empirical  model  is 
conveniently  used  to  determine  conditions  resulting  in  maximum  or  minimum  power 
values.  Furthermore,  udog  plots  of  the  so-called  power  trace,  it  is  possihle  to  detect 
inJiuentia]  error  variances  and  assess  their  impact  on  the  power.  The  power  of  the 


error  variances  is  formulated  in  Section  4.3.2.  In  Section  4.3.3,  the  {tower  value  is 
modeled  udng  generalized  linear  models  tociiniques.  Tb  do  this,  a imrniter  of  designs 
Fue  generated  using  Khuri's  (1996)  method.  Canonical  analy^s  is  performed  in  Sec- 
tion 4.3.4  to  determine  the  nature  of  the  stationary  point  and  the  entire  response 
system.  Optimum  values  of  the  predicted  power  is  obtained  in  Section  4.3.3.  Sec- 
tion 4.3.3  is  devoted  Co  introducing  the  power  trace  plot.  Finally,  concluding  remarks 
are  given  in  Section  4.3.7. 


component  of  the  one-way  random  model  with  unequal 


Consider  the  one-way  random  model 


»u  = ) 


(4.16) 


We  note  that  the  error ' 


‘ equal  within  groups,  that  is,  for  fixed  i, ' 
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ui]i?qu&l  for  different  values  of  i.  Model  (4.16)  is  written  as 

V = /ilv  + (ei'.|l„,)or  + €,  (4-17) 

where  N = U,  Is  a column  vector  of  n,  ones  (i  = 1.2, or  = 

(ouO},"  • , Ok)',  and  y and  e are  the  vectors  of  observations  and  random  errors, 
r^ectively.  The  variancfHurvarlancc  matrix  of  v<  denoted  l^  S,  is  of  the  form 

E = trj  J„,  + ef=,(irf  To,).  (4.18) 

The  usual  F<stalistic  for  testing  the  hypothesis  is 


vrhere  Q and  R are  defined  in  (3.7)  and  (3.8),  respectively.  This  test  is  based  on  the 
assumption  that  the  error  variances,  that  is,  the  cr’'s,  are  equal.  In  this  case,  F has 
under  Wo  the  /’-distribution  with  b - 1 and  W - b degrees  of  freedom. 

Let  us  consider  the  power  of  the  /’-lest  in  (4.19)  when  in  reality  the  error  variances 
ace  not  necesborily  equal.  In  this  case,  the  variaiice-covariance  matrix  C has  the  form 
given  by  formula  (4.18).  The  power  function  is 

P[F  > 0)  = P{y’{Q  - fi,R)y  > 0|ol  5^  0],  (4.20) 

where  A,  = and  is  the  upper  o-percentile  of  the  F distri- 

bution with  b - 1 and  W — b degrees  of  freedom.  Smee  y is  normally  distributed,  the 
quadratic  form  y'{Q  - h«A)y  can  be  cxprtnsed  as 

V’(0-ft,R)V  = |;A,x;,.  (4.21) 

where  Ai.  Ae,  - .• , A,  are  the  distinct  nonaero  eigenvalues  of  (Q  - AoR)E/o*  with 
multiplicities  • • • ,tv,  respectively,  sudi  that  2^,0,  = rank  (Q  - /inR),  and 
the  xj),’s  are  indepeodeut  central  chi-squared  variates  with  or  degrees  of  freedom  (see 
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Johnson  and  Kotz,  1970»  p.  153).  Note  that 


2 = J„,+ 


where  p,  — , i - 1,2,  • 


Fhr  a jiven  k.  the  values  of  and  v,  (<  = 1,2,---  ,r)  in  (4.21)  depend  on  the 
.4NOVA  design  D = {ni,  nj,  - - - , n*}  and  on  the  ratios  /),  (i  = 1,2,--  - , fc).  For- 


mula (4.20)  can  then  be  written  as 


>0 


(4.22) 


Forgiven  k and  jV.  the  power  value  in  (4.22)  requires  the  speciOcotion  of  the  ANOVA 
design  O and  the  ratios  a (ss:  1,2, ,h).  The  right-hand  side  of  (4.22)  can  be  easily 
evaluated  using  the  computer  algorithm  by  Davies  (1930).  The  degree  of  imbalance 
of  the  design  D is  determined  by  a measure  of  imbalance  denoted  by  4,  which  was 
defined  in  (4.3). 

In  the  next  section  we  show  how  to  develop  an  empirical  relationship  between 
the  power,  computed  using  (4.22).  and  N.d.pi.pj, • • • ,p».  The  purpose  of  such 
a relationship  is  to  facilitate  the  study  of  the  combined  elfeels  of  imbaJanen  and 
liecerogeneity  of  the  error  variances  on  the  imwer  of  the  F-test  in  (4.19).  A hey 
ingredient  in  the  development  of  this  relatioosbip  b a method  for  the  generation  of 
an  ANOVA  design  D = {ni,  nj,  ■ ■ - , n*}  for  specified  values  of  N and  We  refer 

(1996).  Note  ciiat  several  ANOVA  designs  can  be  generated  using  the  same  values  of 
.V  and  ^ since  N and  d do  not  determine  the  ANOVA  design  D = {rii.na,  - - - ,ns} 


143 


■1.3.3  Modeling  rhg  pciwpr  of  the  F-tpat 

Suppose  that  values  for  JV.  0,  pi,  - ■■  have  been  specified.  Using  Khiiri'a  (1996) 
method,  several  AN'OVA  designs,  D ss  {ni.nj,'"  ,n*},  can  be  generated  for  Che 
given  values  of  JV  and  d.  For  each  generated  £),  Che  values  of  ni.ns,  • • • ,nt  and 
the  selected  values  of  pi,pe,  ■ ■■  ,ps  are  used  In  formula  (4,22)  to  compute  the  power 
value  which  we  denote  by  rr.  We  can  therefore  regard  a as  a response  vari^le  and 
/V,^,P],p3,  ••  • ,Pa  as  Its  control  variables.  Since  several  ANOVA  designs  D are  gem 
crated  for  eadt  choice  of  N and  d,  a number  of  replications  on  rr  are  available  at 
the  point  (lV,d,Pr,ps,***  ,pv)-  These  replications  will  be  useful  in  the  construction 
of  such  a relationship.  Given  the  nature  of  the  response  rr,  it  would  be  appropri* 
ate  to  model  rr  against  lV,d.pi,P3.‘ " ,pt  using  generalized  linear  models  (GLMs) 
techniques. 

Tb  demonstrate  the  development  of  an  empirical  model  relating  a to  N,  d,  and 
Pi,  ‘ ‘ ‘ >pk,  let  us  cutimder,  as  an  example,  that  b = 6.  Thus,  D = {ni,n2,na,ni,r^}, 
Let  r denote  the  point  [N^  d,  pi , ps,  ■ ■ ■ , ps]^  For  each  selection  of  r In  some  region  of 
interest,  several  determinations  on  C are  made  using  the  corresponding  values  of  N 
and  d in  r.  The  values  of  a are  then  computed  using  Davies'  algorithm  and  formula 
(4.22).  This  process  results  in  several  replicates  on  a st  the  point  r. 

Let  us  now  consider  the  rrsnsformatiou 


GLMs.  For  this  purpose,  we  need  to  specify  a distributioo  for  uj  io  addition  to  a linh 
function,  y,  sucli  that  = g[p(r)].  where  p(r)  denotes  the  mean  of  the  distributioo 
of  at  r.  g is  a strictly  monotone  iocreasing  functioo,  and  { is  a linear  predictor. 
Our  choice  of  g is  the  logarithmic  function.  log(*],  and  that  of  the  distribution  of  ui 
is  the  gamma  distribution.  The  choice  of  this  distributioo  will  be  justified  later.  As 
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for  the  lineer  predictor,  we  consider  a second-degree  model  of  the  form 

e = A4-r'0-fT'flr,  (4.24) 

where  \ and  the  elements  of  d rmd  B arc  unknown  parametera.  Fitting  model  (4.24) 
requirce  usng  on  appropriate  second-degree  response  surface  design.  A convenient 
choice  is  the  Box  and  Behuken  [1960]  design  for  seven  variables  with  four  center-point 
repllcatfs.  This  is  an  economical  3-lcvel  Incomplete  factorial  design.  The  selected 
three  levels  for  cadi  of  the  elements  of  r = (JV.d.pi.pj,  • ,ps)  are;  N — 25,50,75; 
4 = 0,35.0.65,0.95;  p,  = 0.1, 0.5, 0.9,  i ° 1, 2,-'  ,5.  The  actual  settings,  in  coded 
form,  of  the  Bos-Behnken  design  are  given  in  Table  4.12.  Note  that  except  for  i'f,  the 
values  of  all  the  remainiug  elements  of  r fall  inside  the  unit  interval  (0, 1|.  To  reduce 
the  range  of  N,  we  consider  a logarithmic  transformation  of  N.  Xi  = iog(A').  This 
way,  the  range  of  Xi  for  the  selected  values  of  A'  is  about  1,1,  which  is  comparable  to 
the  ranges  of  4 and  ft  [i  = 1,2,--  - , 5).  Let  ij  = 4,  ij  = hi,  is  = pj,  *s  = tb,  = 
ft.  x^  = Pt,  and  « = {i|,i3,  - ,17}'.  Model  (4.24)  can  then  be  expressed  in  terms 

( = ai  + x'B  + x'Bx,  (4.25) 

where  and  the  elements  of  0 and  B are  unknown  parameters. 

For  the  chosen  levels  of  X{  (i  s - ,7),  the  region  of  interest  for  the  fitting  of 
model  (4.25)  is 

C = (x|  log(25)  < i|  < log(75).  0.3  < n < 1-0,  0.1  < i,  < 0.95.  i = 3,  • • • , 7}(4.26) 
Using  the  logarithmic  link  function  in  coraldnation  with  rnodel  (4,25),  we  get 

loglA(*)l  = &,  + x'0*  x'Bx.  (4.27) 

where  p{x)  is  the  mean  of  the  distribution  of  u>  in  (423)  at  the  point  x.  It  should 
be  noted  here  that  the  logarithmic  link  was  preferred  over  the  reciprocal  link,  which 
is  the  canonical  link  Function  for  the  gamma  distribution.  The  reason  for  this  will  be 
explained  later. 


Table  4.12.  BcK-BehnIcen  design*  and  the  number  of  ANOVA  designs  genernted 


d Pi  P3  hs  Ai  Ps  Number  of  designs  generated** 
0 0 0 ±1  ±1  d:l  0 6 for  each  of  S design  points 


fl  (I  ±1  ±1  0 0 ±1  6 for  each  of  S design  points 

±1  0 ±1  0 ±1  0 0 a 

0 ±1  ±1  0 0 ±1  0 b 

DDOOOOO  6 for  each  of  4 center  points 


(a)  12  designs  for  iV  = 1;  6 designs  for  Af  = -1 

(b)  12  designs  for  ^ = 1;  7 designs  for  d = -1 

(c)  20  designs  for  {S,  d)  — (l.l)i  14  designs  for  (Af,d)  = (1,  *-l)i  4 designs  for  each 

of  (N,d)  = (-M).(-1.-I) 


The  coded  settinp,  -1,0,1.  denote  the  low,  middle,  and  high  levels, 
Tlie  total  luimber  of  ANOVA  designs  generated  is  bOO. 
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At  c&ch  point  of  the  Box-Bclmkon  design  in  l^blo  4.12,  several  ANOVA  designs  are 
generated  usig  Khuri's  (1990)  metiiod.  It  gives  severai  replications  on  s,  and  hence  on 
u,  at  each  point  of  tile  Box-Behnken  dedgn.  The  number  of  generated  determinations 
of  iD  = (ni,nj,nj,nj,nj}  corresponding  to  each  point  of  the  Box-Belinken  design  is 
also  given  in  Tible  4-12.  Some  of  these  determinations  along  with  the  ectual  d-volue 
for  each  generated  design  are  shown  in  Table  4.13.  The  corresponding  values  of  a, 
based  on  the  fonrailns  (4.22)  using  Davies'  algorithm,  are  also  given  in  thble  4.13. 
Let  Uf  and  denote  the  sample  mean  and  sample  variance  of  the  w.replicatos  at  x 
corresponding  to  a point  x in  the  Box-Behiiken  design.  The  task  now  ie  to  determine 
if  there  is  a strong  linear  relationship  between  and  or  between  iSt  and  s^.  Fitting 
a sample  linear  regression  model  with  no  intercept  to  the  (C„  s‘)  data  results  in  s 
coefficient  of  determination  fP  = 0.80.  Fitting  the  same  model  to  the  (Wj.J,)  data 
gives  re  = 0.97.  This  indicates  a strong  linear  roiationsliip  between  the  mean  of  the 
distribution  of  w and  its  standard  deviation,  wliicli  justiffes  our  choice  of  the  gamma 
distribution  for  tv. 

Fitting  model  (4.27)  to  the  u-data  that  were  generated  ot  the  points  of  the  Box- 
Behnlten  design  produces  s predicted  value  for  tv,  and  hence  fortr  [see  formula  (4.23)]. 
We  llien  have 

»(*)  5 — 1 — , (4.281 

1 + expos, + *'13  + ST’S*)  ' 

where  w(a)  denotes  the  predicted  power  value  at  a point  z in  the  region  C defiued 
in  (4.26),  and  and  tlie  elements  of  Q and  B are  the  maximum  likelihood  atimates 
corresponding  to  ft,  0,  and  B,  respectively.  The  GENMOD  procedure  in  SAS  (1897) 
was  used  u>  lit  model  (4.27). 

Let  the  i*^  element  of  ^ and  the  (t,;)“  element  ofB  in  (4.28)  be  denoted  by  ft  and 
ft,,  respectively  i = 1, 2,*  • • ,7;  j * 1,2,  - • • ,7.  The  tests  concerning  ft, ft,, ft,, ft,, 
and  ft,  were  found  to  be  highly  insignificant.  These  terms  were  therefore  dropped 


from  the  model.  Model  (4.27)  was  aubseqiiently  refitted  using  the  remaining  pa- 
rameteni.  Estimates  of  these  paramctore  are  given  in  Table  4.14.  We  note  that  all 
pararoeiere  are  significantly  diflerent  horn  sera.  With  theae  parameter  etimates,  the 
estimated  power  values  at  all  points  of  the  Box-BehnJien  de^  were  computed  using 
model  (4.28).  The  results  areshown  m Table  4.15.  Since  there  arc  several  replications 
on  7T  at  each  design  pomt,  the  maximum  dilference,  as  well  as  the  averse  difference, 
between  the  replicated  values  and  the  corresponding  estimated  value  of  s (at  the 
same  location)  Is  computed  to  check  the  adequacy  of  fit  of  the  model.  These  values 
are  also  pven  lo  Table  4.15.  It  can  be  noted  that  the  maximum  difference  raoges 
from  -0.1077  to  0.1296,  except  at  one  point  where  it  is  equal  to  -0.2412.  The  values 
of  average  dilference  range  from  -0.0748  to  0.0916,  except  at  one  point  where  it  is 
equal  to  -0.2073.  These  values  along  with  those  of  the  scaled  deviance  and  scaled 
chi-squared  statistic  in  Thbic  4.14  indicate  that  mode]  (4.2S)  provides  a good  St  to 
the  data. 

It  should  be  noted  here  that  using  the  canonical  link  function  for  the  gamma  dis- 
tribution, namely,  the  reciprocal  link  function  instead  of  the  logarithmic  link  function, 
produced  Infeasible  rosulla:  ai  some  points  of  the  Box-Belinken  demgn,  the  estimated 
value  of  s did  not  fall  inside  the  inlervnl  (0,  Ij. 

Note  that  the  ranges  of  interest  for  ij  (Le-,  for  ^)  and  Jj,  j = 3,4,  • • • ,7  (i.e.,  for 
A,  i = 1,2,  - ,5),  in  (4-26)  cover  wide  range  of  the  parameter  apace  ofi.’s,  namely, 
0 5 S 1- 1 = 2,3,  • - - , 7.  Tb  check  the  prediction  capability  of  the  model  (4.28)  at 
points  outside  the  region  C - especinlly  for  the  variable  xi  (le-,  for  jV)  - we  generate 
designs  with  selected  values  of  x.  These  demgns  along  with  the  values  of  x are  listed 
in  Table  4.16.  Note  that  iV  = 15  and  N = 100  are  the  values  whidi  fall  outside 
the  region  C in  (4.26).  Also,  the  value  of  0 = 1,  whidi  is  for  the  balanced  design 
case,  was  not  used  to  build  the  model  (4.28).  Hence,  its  prediction  capability  of  the 
model  (4.26)  whea  d = 1 is  also  of  interest. 


Tabic  4.13:  Some  generaied  ANOVA  dcagns  and  the  cor- 
responding  Cnie  v^ues  of  ir  denotes  the  actual  value 
of  ip  for  a ^generated  design) 


and  gocdness-of-fit 


ISO 
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Table  4.15.  — continued.  The  estimated  power,  lr(x),  and  Its  maximum  and  average 
differences  with  xj,  the  i"  power  value  at  a given  point  of  the  Box-Behnlien  design. 


For  each  design  in  Table  4.16,  the  exact  power  is  calculated  using  Davies’  algo- 
rithm. The  predicted  power  is  also  obtained  uuug  the  model  (4.25).  These  valus 
along  with  its  differences  are  listed  In  Table  4.17.  f>om  this  t^e,  we  can  see  that 
the  prediction  capability  of  the  model  (4.25)  is  still  good,  even  at  points  uotaide  the 
region  C in  (4.26). 

4,3.4  Canonical  analysis  of  model  14.281 

Model  (4.25)  with  paranieter  estimates  in  Table  4.14  can  be  used  to  determine 
conditions  on  z = (ii.zj,  - - - ,17)'  that  result  In  large  values  of  if.  We  recall  that 
ii  = log(lV),  ij  = ij  = Pi,  Is  = P2,  *6  = ps,  ij  = p,.  x,  = pj.  Let  0,  and 
^ denote  Che  lower  and  upper  bounds,  respectively,  on  i,*  as  described  in  the  region 
C in  (4.26),  i = 1,2,  • • • , 7,  It  is  convenient  here  10  code  the  i|’s  using  the  variables 
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Table  4.16.  Generated  designs  Tot  one-way  model  with  1:  = 5 to  check  the  prediction 
capability  of  model  (4.3S).  denotes  the  actual  value  of  d for  a d-generated  design 


Table  4.17.  Prediction  capability  of  model  (4.28)  for  designs  in  Table  4.16 


i = Gi  + d (4.30) 

where  * = (r„  r,,  • • • .s,)',G  = ®L,  (j^),  end  the  element  of  d is  d,  = Jit, 

< = T»  + *'7  + i'r*,  (4.31) 


where  To  = A - d'G-‘/3  + d'G-‘BG-'d,  7 = G-’/3  - 2G-'BG-‘d,  itnd  T = 
G-'BG-'.  Theestimnied  value  of  {is 

j = iu  + *'7  + »'ft.  (4.32) 


and  B.  FVom  (4.28),  the  predicted  value  oftr  at  a point  » in  the  re^on  C„  where 


C.  = {*|-l<e,<l,  : = l,2,.--7}.  (4.33) 


is  given  by 


T7K)- 


• off; 
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ortLogun&l  matrix  that  dfagonaliaes  t,  that  ja,  f - PAF*.  where  A is  the  diagODal 
matrix  of  eixenvaluea  of  F (sec  Kliari  and  CoroeU,  1996,  Section  5.6.1).  Note  that 
the  cotimms  of  P are  orthonormai  eigenvectors  of  F coirespondiog  to  Kg,ics,  • • • , K?. 
Using  the  parameter  estimates  in  Table  d.l4,  we  have 

& = U-4817 

So  = (-2.638.3.082,-9.302,-9.408,-6.081,-13.140,5.094)' 


A = diag(1.370.0.612.0.421, -0.400,0.206, -0.041, 0.020).  (4.30) 


We  note  that  2g  falis  outside  the  region  C,  io  (4J3).  Aiso,  since  the  eigenvalues 
in  (4.36)  have  variabie  agns,  sg  is  a saddle  point. 

The  main  advantage  of  the  canonical  equation  (4.35)  is  the  ability  to  study  the 
behavior  of  ir  along  the  principal  IV, -axes,  that  is,  those  axes  that  are  determined  by 
the  roves  of  P".  This  is  done  as  follows;  since  a,  and  >t«  are  negative  [see  (4.36)|,  a 
move  away  from  in  along  the  HVaxis  or  the  Hi-axls  in  either  direction  resulii  in  n 
decrease  in  { nnti  hence  an  increase  in  ir  since  ir  = ^ ^ . The  increase  in  i is  more 
sizable  along  the  lV,.axis  since  |X||  is  larger  than  |xa|.  On  the  other  hand,  a move 
away  from  «o  along  each  of  the  (t'faxis,  i - 1, 2,3,5, 7,  in  either  direction  results 
in  an  increase  in  i and  hence  a decrease  in  i.  The  decrease  in  ir  is  most  apparent 


ftlong  ihe  H'l-axis  since  Ki  is  larger,  in  terms  uf  absolute  value,  than  alJ  the  other 
eigenvaJues- 

Since  the  stationary  point  zo  is  a saddle  point  and  foils  outside  the  region  C| 
in  (4.33),  we  should  detennine  the  absolute  optima  (maximum  and  minimum)  of  x 
within  C„  or  equit^ently,  the  region  C in  (426),  This  issue  wili  be  addressed  in  the 


Let  118  write  the  linear  predictor  in  (4.25)  as 

£ = (4.37) 

where  4 is  a vector  of  pararaetere  and  /'(at)  is  a vector  of  powers  and  cross  products 
of  powers  of  the  elements  of  e up  to  degree  2 such  that  /'(®)4  = Bi  + x‘0*  z‘Bx. 
Dy  evaluating  x at  each  point  of  the  Box-Behnken  desip  (mentioned  earUer  b Sec> 
tion  4.3.3)  wu  get  tlie  model 

« = X4,  (4.38) 

where  X = |/(*i)  : /(wj)  ; - : /(«:„)]',  n = S,m(  = 500  with  m,  being  the 

number  of  replications  at  the  desip  pomt,  i s 1,2,  • • • ,60,  and  ( is  the  n x 1 
vector  of  corresponding  values  of  J in  (4.37).  The  predicted  value  of  f at  a point  x 
ill  the  region  C in  (4,36),  denoted  by  4(s:),  is  ((a)  — /'(st)4,  where  4 is  such  that 
/'(»)4  = ,^  + x'0  + x'Bx.  The  variance  of  4(tc)  is  approumately  given  by 

Var[{(®)|  = /'(x)(X'A-'X)-7(*).  (4.39) 

where  A is  a diagonal  matrix  whose  t^  diagonal  element  is  ^[t;’(ii,)]^,  i = 1,  • • • ,n, 
where  p(ii,]  = log(fii)  with  p,  being  the  mean  of  the  distribution  of  u in  (4.23)  at  the 
desip  point,  and  p'(p,}  is  the  6rst  derivative  of  g evaluated  at  /a.  Furthermore, 
by  our  choice  of  the  gamma  distribution  for  tv,  the  variance  of  this  distribution  at 
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ilic  point  i(i  ^ven  by  where  v is  a scale  parameter.  An  estimate  A of  A is 
obtained  via  an  ilerativeiy  weighted  least-squares  procedure.  Hence,  an  estimate  of 
the  variance  of  i(*)  = is  approximately  given  by 


The  function  9(tc)  is  a scaled  version  of  x(s:)  and  is  the  one  Co  be  optimized  over 
the  region  C.  This  is  accomplished  numerically  through  a grid  search  over  C.  It  was 
found  Chat  q(x)  attained  its  minimum  value,  namely  5.174,  at  Che  point  x„,„  whose 
coordinates  are  xi  = log(/V),  with  N = 25.  jj  = ^ = 0.3,  Zs  = pi  = 0.1,  xu  s 
P2  ~ 0.95,  xb  = Pa  = 0.95,  Xs  = Pt  = 0.95,  and  x?  s*  = 0-95.  The  com«ponding 
value  of  i(*)  is  0.004  with  a standard  error  approximately  equal  to  0.00078.  The 
maximum  value  of  q(x)  over  C is  15184.31.  and  is  attained  at  the  point  whose 
coordinates  are  Xi  = log(^),  with  IV  = 75,  is  = d = 1-0,  Xs  = pj  = 0.95,  xt  w p,  w 
0,95,  xj  = Pa  = 0.95,  is  = pi  = 0.95,  and  17  = ps  = 0.95.  The  corresponding  value 
of  i(x]  is  0.99956  with  a standard  error  approximately  equal  to  0.000066.  We  note 
that  both  optima  occur  on  the  boundary  of  the  region  C.  We  also  note  that  in  Xm<n, 
the  values  of  jV.d,  and  pi  are  equal  to  their  lower  bounds  in  C,  whereas  the  values  of 
Pg.fl7.Pi,  and  fig  are  equal  Ui  their  respective  upper  bounds.  On  the  other  hand,  tiie 
values  of  W.d,  and  pi  in  Xmw  coincide  with  their  upper  bounds,  and  pt.ps.pi,  and 
Ps  maiiitaiu  the  same  values  as  in  Xm.,,.  Thccofore  the  maximum  power  occurs  when 
the  errois  are  homoscedastic  (recall  that  p,  — the  data  set  is  balanced  (since 
d = 1),  and  Che  value  of  jV  is  largest  within  the  region  C in  (4.26). 


;/'(x)[X‘A-'x|-'/(x).  (4.40) 


Let  us  now  divide  x(x)  by  its  standard  deviation  to  get 


{Var(x(x)i}» 


(4.41) 


]n  this  section  we  show  how  to  monitor  the  increase  in  the  power  values  as  we 
move  from  *„„,  to  This  is  accomplished  by  uring  plots  of  the  so-called  power 
tracer  which  we  now  describe. 

U&Dg  the  coding  convention  of  formula  (4.29),  Smm  and  can  be  written  as 
*nin  = 1)',  tma  = (1, 1. 1, 1, 1, 1, 1)',  respectively.  Let  us  trans- 

late the  2-coordinaiS5  system  so  that  its  origin  coincides  with  s„,,„.  This  is  accom- 
plished by  the  transformation  u - a - x„,„,  where  a = (zi.sj,*  • • , s?)'.  FVom  (4  JO), 
the  relationship  between  x and  a is 

x = G-’(u-d-ha„„).  (4,42) 

The  region  C,  in  (4.S3)  can  then  be  expressed  ns 

C,  = {u|Q<it.<2,  i = l,2,3;  -2  < ik  < 0,  i = 4,5,6,7}, 

(4.43) 

where  u,  is  the  i**  element  ofu,f  = 1,2,--*  ,7.  The  points  a„,„  and  Bmoj  corrtspond 
to  u„,,„  s (0, 0, (1, 0, 0, 0, 0)’, w (2.2,2,0,0,Q,0)‘,  respectively.  In  addition,  the 
scaled  predicted  power  function  «(«)  in  (4.41)  can  be  written  as p(u),  where  by  (4.42), 


p(u)  = 7lG-'(u-d 


In  order  to  study  the  increase  in  power,  that  is,  in  Che  values  of  ^u),  as 
from  UnQn  to  UjTiMs  within  the  region  C,.,  we  propose  to  hod  maximum  values  of  p(u) 
on  the  surfaces  of  concentric  hypempheres  of  varying  radii  centered  at  Unm*  This 
is  analogous  to  the  method  of  ridge  analysis  commonly  used  in  responses  surface 
methodology  (see,  for  exsmplc,  Khuri  and  Cornell,  1996,  Section  3.7).  In  our  case. 
maximisatioD  is  confined  to  the  purtitms  of  the  hyperspheres  contained  inside  C^.  This 
process  is  facilitated  by  expressing  ui.us,  • - - ,Ut  in  terms  of  the  spherical  coordinates. 


tl;  = $5in^C05^ 

U3  = dsinti'i  slnti's  C091/9 

ti4  s Aeini^i ^nV'ssui^cos^t 

115  = AsiaiAi^i^sinV^an^iOKi/^ 

Us  = stdostisjns£)38iTn^einV's6in  tt'scost^^ 

U7  = AHnV'|8in^6iDi/^sui  ^sSbsE^^i^, 

where  0 < • < 2V3,  0 < iSr  < if,  i = 1,2, 3, 4,5;  0 < «^  < 2i.  Here,  2^3  is  the 
disumce  of  u^m  from  u,n,n. 

For  a given  radius  s,  points  are  generated  on  the  surface  of  the  hypersphere 
T,  — ,{u|u'u  = by  selecting  ,da  at  random  using  independent  uni- 

form distributions.  Points  that  do  not  fall  inside  are  discarded.  The  remaining 
points  ate  used  to  ratimate  the  maximum  of  p(u)  on  T,.  Let  such  a maximum  be 
denoted  by  ptr^{s).  The  point  at  whitii  this  maximum  is  attained  is  denoted  by 
= [u5'*,t4’*, •'•  Tfif  neitt  step  is  Ui  plot  Pm„(s)  as  web  as  ti[‘*  against 

s,  i = 1,2,  • • • , 7.  Tiiese  plots,  which  we  refer  to  as  power  trace  plots,  are  shown  in 
Figure  4.7.  Note  that  Figure  4.7(a)  gives  a plot  of  the  scaled  predicted  power.  q{x) 
in  (4.41),  against  s,  and  Figure  4.7(b)  is  for  the  corresponding  predicted  power,  w(x) 
in  (4i8),  against  s.  Also,  Figure  4.7(c)  gives  plots  of  ij'*  against  s,  where  j!*’  Is  the 
element  of  zE^  - ulnis  + Zntfl.  The  reason  far  using  a**^  instead  of  is  due  to 
the  fact  that  -I  < Sj'*  < 1 for  i = 1, 2,*  • • ,7,  which  faciiitates  the  supprimposilion 
of  the  seven  plots  of  2)*'.  These  plots  correspond  to  log(N),0.pi,pa,p|,p7,  and  flv. 


rt'spectively. 


161 


From  Figure  4. 7(b,c)  it  can  be  easily  seen  that  pi,ps,pi,  and  /)»  exert  little  eScct 
on  the  power.  On  the  other  hand,  N,i,  andpi,  in  particular,  contribute  sgnificantly 
to  the  increase  in  power-  Thue,  the  power  is  more  sensitive  to  changes  in  N.d,  aiidpi- 
Note  that  for  s < 2.4,  p\  is  moat  inhuentiai,  roliosod  by  ^ aud  N.  Hence,  the  power  is 
most  sensitive  with  respect  to  pi.  We  mx)  recall  that  pi  is  associated  with  the  largest 
group  sice  ni  since  the  n,  values  in  Tabie  4.13  were  arranged  in  descending  order. 
Figure  4.7(c)  shows  that  the  values  of  pi  ercarmillet  than  those  of  p,  fori  = 2, 3, 4, 5, 
hence  af  is  larger  than  trf  far  i = 2, 3, 4, 6.  We  conclude  that  aligning  the  largesrt 
number  of  experimental  unita  to  the  moat  variable  group  raaltes  its  variance  most 
Influential  with  regard  to  the  power  of  the  F-test.  Thus  assigning  the  largest  sample 
mo  to  the  group  having  largest  variance  can  produce  high  power  values  provided  that 
the  iargesi  variance  Is  not  too  large,  noting  that  if  the  largest  variance  is  small,  the 
other  variances  arc  also  small. 

4.3.7  Concludine  remarks 

The  modeling  of  the  power  of  the  F-test  against  fV,  pi , pj,  • • • , p»  enables  one 
to  clearly  assess  the  effects  of  these  parameters  on  the  power.  Response  surface 
techniques  can  be  couvertiently  used  for  this  purpose.  In  particular,  the  canonical 
analysis  in  Section  4.3.4  is  utilized  to  determine  the  nature  of  the  statiooaiy  point 
of  the  fitted  surface,  and  to  identify  principal  axes  along  with  the  power  increases  or 
decreases.  Funhermore,  plots  of  the  power  trecc  in  Section  4.3.6  provide  a graphical 
depiction  of  the  increase  in  the  power.  These  plots  are  also  useful  in  detecting  those 
error  variances  that  are  influential  contrtbutoiii  to  the  increeae  in  the  power. 
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4 j The  Covertte  ProtabililY  of  ft  ConfidOTce  InttrvAl 
4.4.1  Introduction 

Confidence  intervals  on  variance  components  are  of  interest  in  any  study  which 
deals  with  random  or  mixed  effects  models.  Many  papers  have  been  written  on 
the  subject  of  confidence  btervals  for  variance  components,  functions  of  variance 
components,  and  simultaneous  confidence  intervals  for  a variety  of  models.  For  a 
wide  coverage  of  this  subject  area,  we  refer  to  Burdick  and  Graybill  (1992). 

An  exact  confidence  interval  on  a function  of  variance  components  does  not  always 
exist.  This  is  true  even  for  the  one-way  random  model  with  two  variance  components, 

and  of.  t^Tien  the  dataset  is  balaoced,  there  are  exact  confidence  intervals  on  a*, 
the  ratio  and  the  intraclass  correlation  coefficient,  p — (see  Burdick  and 
Graybill,  1992,  Ch.  4).  However,  there  does  not  exist  an  exact  confidence  Interval  on 
o|  as  Thomas  and  Hultquist  (197S,  p.  3S3)  Indicated.  Several  studies  were  made  to 
construct  an  approximate  confidence  interval  on  (see,  for  example,  Bartlett,  1953; 
Bross.  1950;  Bulnier,  1957;  Cochran,  1951;  Howe,  1974;  Moriguti,  1954;  Satterthwaite. 
1941;  Tiikey,  1951;  Welch,  1956;  Williams,  1962,  to  name  just  a few).  Boardman 
(1974)  compared  those  apprcrdmalc  procedures  by  using  Monte  Carlo  simulation  to 
determine  the  actual  probability  coverage  for  each  procedure.  His  study  showed  that 
two  procedures  provided  a coverage  probability  close  to  the  nominal  level.  One  of 
these  procedures  was  independently  develoited  by  'Bikey  (1951)  and  WilUaros  (1962). 
The  other  procedure  was  independently  developed  by  Moriguti  (1954)  and  Bulmer 
(1957). 

For  the  unbalanced  one-way  random  model,  only  of  has  an  exact  confidence  inter- 
val. For  ocher  parameters  of  interest,  several  approximate  confidence  intervals  have 
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been  proposed.  Foremmplc,  foro|,  the  work  of  Thomas  and  Hultquist  (1978).  Bur- 
dick and  Eickmim  (1886),  and  El-Baasiouni  (199-1)  are  noted.  Greybill  and  Wang 
(1980)  considered  a balanced  ff-factor  nested  random  model.  They  proposed  a gen- 
eral procedure  for  constructing  a confidence  interval  on  a linear  comblnntlon  of  vari- 
ance components  based  on  K Independently  distributed  chi-squared  random  variables. 
Burdick  and  Craylnll  (198d)  proposed  a modification  of  the  Graybill-Wang  procedure 
that  considered  the  unbatancedness  of  the  design.  Their  modification  was  for  A'  = 2 
whidi  is  for  the  unbalanced  one-way  random  model.  Burdick,  Maqsood  and  Gray- 
bill  (1986),  Donner  and  Wells  (1986),  and  Mian,  Shoukri,  and  TVacy  (1989)  obtained 
confidence  intervals  on  the  ratio  of  variance  components,  and  the  intraclass  corre- 
lation coefficient.  For  models  other  than  the  one-way  random  model,  Kaaempour 
and  Grayblll  (1991)  considered  the  unbalanced  two-way  random  cross-classification 
with  interaction  modeh  Gcaybill  and  Wang  (1979)  and  Wang  and  Greybill  (1981)  fo- 
cused on  the  balanced  two-fold  nested  random  model;  Sen,  Gcaybill,  and  Ting  (1992), 
Hemandes,  Burdick  and  Birch  (1992),  Hernandez  and  Burdick  (1993)  considered  the 
unbalanced  two-fold  nested  random  model;  and  Harville  and  Fenech  (1985)  dealt  with 
an  nobalanced  mixed  linear  model.  For  a survey  on  confidence  intervals  up  to  1984, 
see  Khuri  and  Sahai  (1985). 

Even  though  many  papem  were  written  on  confidence  intervals  for  variance  com- 
ponents, there  is  no  study  devoted  mainly  to  the  effect  of  Imbalance  on  the  coverage 
probability  of  confidence  intervals.  The  coverage  probabiUty  of  a ~ri.iT,  confidence 
interval  depends  not  only  on  the  variance  components,  but  also  on  the  dmign  used. 
Researchers  have  used  particular  designs  ranging  from  extremely  unbalanced  to  neatly 
balanced  in  order  to  evaluate  the  performance  of  their  proposed  methods. 


Ill  Seciioji  iA,  we  shall  invesiigate  the  p|?ccl8  of  imbalance  and  variance  compo< 
nenteon  the  coverage  probabilities  of  several  confidence  intervals.  Fbt  the  unbalanced 
one-way  random  model,  four  approximate  methods  of  confidence  Interval  on  are 
considered.  These  are  given  in  Section  4.4.2.  In  Section  4.4J,  we  ptopoee  a proce- 
dure for  modeluig  the  coverage  probability  associated  with  each  of  these  methods  in 
terms  of  factors  that  infiuenee  the  coverage  probability.  Section  4.4.4  is  devoted  to 
com]Sorlng  these  methods.  Concluding  remarks  are  given  in  Section  4.4.fi. 

4.4.2  Confidence  intervals  on  o!  for  the  Qiie-wav  model 


Consider  the  one-way  unbalanced  random  model  (3.1)  and  its  vector  form  in  (3.2). 
The  sums  of  squares  for  the  treatment  effect,  and  the  error  term.  y'Ry,  where 
Q and  R are  defiiied  in  (3.7)  and  (3.S),  respectively,  are  independent.  When  the 


arc  not  equal,  y'Qy  has  a scaled  cht-squared  distribution  if  and  only  if  = 0.  If  It 
is  known  that  is  close  to  xero,  it  may  l>e  appropriate  to  treat  i/'Qv  approximately 
as  a scaled  chi-squared  variable.  This  motivates  the  use  of  an  apiiroxiraate  confidence 


Modified  Large  Sample  (MLS)  Confidence  Interval 

Burdick  and  GraybiU  (1992,  p.  fil,  p.  70)  formulated  a conSdeoce  interval  based  on 
the  above-mentioned  strategy.  The  resulting  approximate  two-sided  (1  — o)100% 
confidence  interval  on  is 


ffi/So  - AfSe  - vTT  fifS.  - AfSg  + vAtil 


1S5 


MS,  = 

.WSs  = ij^y'Ry. 

\\  = C?.WS2  + H|A#5|  + G,jAfS..W5e, 

Vc  = H^MSi*ClMSI:  + H„MS,MSs. 

* /|(k-l,oo)' 

“ F|.,(*-l, «.)■*' 

. |f,(*  - 1,  Af  - t)  - 1]»  - |G,f,(t  -\.N-  k)]'  - fl? 

■'  F,{k-\.N-k) 

[1  - - 1,  A-  - *)I»  - 1,  A - i)l=  - d 

Since  > 0,  tbe  lower  bound  Is  set  to  zero  if  < F^{k  — 1,R  - k).  Also,  tbe 
upper  bound  is  defined  to  be  scro  if  < /i-}(fc  - l.N  - *).  Here, 
denotes  the  upper  IDOd*^  percuntile  of  an  F-distribution  with  vi  and  degrees  of 
freedom.  We  shall  refer  to  formula  (4.4d)  as  tbe  MLS  confidence  interval  on 

Thumas-Huliquisi’s  jTH]  Confidence  Interval 


The  unweighted  sum  of  squares  associated  with  the  treatment  means  io  an  unbalanced 
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“heteft.  = = S Ef=i!'.- 


li  is  itolsd  chat 


(4.45) 


is  the  sample  t’arianee  for  the  treatraent  means.  Also  note  that  .5(^)  — 


1964.  p.  152).  Thomas  and  Hnitquist  (1976)  showed  that  the  moment  generating 


squared  random  variable  as  all  rq  values  approach  a constant,  or  as  all  n,  -»  oo,  or 
as  q ^ ^ » 00.  They  also  provided  simulation  results  that  IV  is  wetl  approximated 
by  a chi-squared  random  variable  when  rj  > 0.26  and  the  design  is  not  extremely 
nnltalanced.  Using  IV,  they  derived  forinniae  for  constructing  conhdence  intervals 
on  o^,  |7,  and  ■ More  speciKcally,  on  approximate  two-fflded  (1  — a)100% 
txmfidence  interval  on  ffl  is  obtained  by  replacing  MSa  with  n/,Sg  and  n with  ns  in 
the  WiiliBins-Ibhey's  balanced  model  formula  (Williams.  1962;  Ibkey.  1961).  The 
reaultiug  interval  is 

[tt,5;  - AfSeFfil!  - 1,  A'  - *)  ruSj  - A/SsFi-j  (*  - 1,  A"  - *)  1 

[ n,F5(k-l,oo)  ' n,f|.j(i-l,eo)  J' 

By  replacing  Fj(i-l,cx>)  with  j^y|(t-l)  and  F|_j(i-l,oo)  with  jijj:?.  }(*-!). 
formula  (4.46)  is  equivalent  to  the  formula  given  by  Thomas  and  Huitquist  (1978. 
p 588).  The  lower  and  upper  bounds  of  (4.46)  are  defined  to  be  zero  if  < 
Fj(i-1,  A -k)  and  ^ < F|_j(k-  l,A-k),  respectively,  since  cj  > 0.  Thomas 
and  Hultquist  mentioned  that,  bisituatious  where  au  approximation  of  the  IV-random 
variable  works  well,  this  formula  has  a good  coverage  probability.  We  shall  refer  to 
formula  (4.46)  as  the  TH  confidence  interval. 


where  n#  - 


is  the  harmonic  mean  of  the  n,  values  (see  Burdick  and  GraybiJI, 


function  of  Iheir  IV-randoin  variable,  which  is  W = approaches  that  of  chi* 
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Btirdick-Eickinaii'B  (BE)  Confidence  [nUirval 

An  eppcoximelc  imerva!  is  said  to  be  conservative  If  its  coverage  probability  is  larger 
ihon  the  nominal  value,  and  liberal  if  it  is  smaller  than  the  nominal  value.  In  ex* 
tremely  unbalanced  designs  where  < 0.25,  the  Thomaa-Hultquist  confidence  in- 
terval on  o1  can  be  liberal  (Burdick  and  Graybill,  1992,  p.  70).  The  approximate 
(1  - q)100%  confidence  interval  on  based  on  Thomas-HulUprist'a  approx* 

itnation  can  also  be  liberal,  and  Burdick,  Maqsood  and  Graybill's  (1966)  confidence 
interval  on  ^ is  conservative  since  It  guarantees  a coverage  probability  at  least  as 
large  as  I- a.  Combining  these  two  intervals  using  Willituos'  (1962)  strategy  for  con- 
structing a confidence  interval  on  in  the  balanced  design,  Burdick  and  Eickmau 
(1986)  Introduced  a confidence  interval  on  for  the  unbalanced  one-way  random 

[ ] (4  47) 

[f,(lr-l,oo)(l-t-Tt.i„)'  f,.,(*-l,»)(H-ns£/«)J’  ' ' 

, , S} ^ 

" F,(k-I,jV-*)MJi-  :!(„• 

S 3 

f|.j(k  - l,iV-k)M5e  ri(s)’ 
n^i]  = min  (ni , ri3,  • • - . ne). 
n(*3  = max(ni,n2,' • • ,nk), 

and  rir,  and  S|  are  defined  as  in  (4.45),  Again,  the  lower  and  upper  bounds  of  (4.47) 
as  the  BE  confidence  interval  on 
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Modified  Harmonic  Mean  (MHM)  ConHdence  Interval 

Consider  the  matrix  K = ding(;^.  The  matrix  j Js)  has  fc  — 1 

positive  eigeovalues,  namely,  Ai,  A;,  - • , Ak-|.  The  average  of  these  eigenvalues  is  Che 
same  as  the  reciprocal  of  the  harmonic  mean  of  the  m-values,  i.e.,  A = A,  = 

j Khuri  (1999)  proposed  replacing  114  In  the  TH  conlidence  inter- 
val (•1.46)  with  the  value  n„,  - where  A(i)  and  A(k.i)  are  the  minimum  and 

maximum  of  A,*a.  respectively.  The  resulting  interval  on  is 


and  <F|.f(h-l,A^-h),  respectively.  This  choice  of  On  is  supposed  to  improve 
the  approximation  of  IV  in  Thomas  and  Kultquist  (1976)  with  the  chi*squared  variate 
\l~  I'  refer  to  formula  (4.49)  as  the  MHM  conhdence  Interval. 

4.4.3  Modeling  the  covtrMf  nrilUWIliY 

The  coverage  probability  of  a conddence  interval  on  with  a nominal  value  1 - a 
depends  not  only  on  the  variance  components,  o*  and  a’,  bat  also  on  the  design, 
V,  that  is  the  set  of  n,-vaiues  for  model  (3.1).  Therefore,  it  is  of  interest  to  study 
the  elfecrs  of  these  factors  on  the  coverage  probability  of  a confideoce  Intervzd.  The 
effects  of  variance  componenta  on  the  coverege  probability  have  been  studied  using 
a variety  of  methods  by  several  rceeorchcre  ae  wc  have  seen  in  Section  4.4.1.  In  all 
of  these  studies,  only  a few  designs  that  were  moderately  or  extremely  unbalanced 
were  used.  These  designs  did  not  cover  a wide  range  of  unbalanredness.  Moreover, 


it  was  hard  to  detect  the  rumbmed  effects  of  variance  componeots  and  design  on  the 


Table  4.18.  Number  of  desisns  generated  for  each  combination  oF  N,  0.  and  p 


As  was  seeu  earlier  iti  Chapter  4,  a design  D for  modei  (3.1)  can  he  generated  by 
specifying  the  total  number  of  obsertatlon,  /V,  and  the  degree  of  Imbalance,  tp,  for 
a fixed  number  of  treatment,  k,  usiog  Khuri's  (1996)  method.  Such  a design  is  not 
uniquely  determined.  In  this  section,  we  provide  an  empirical  model  to  represent  the 
txiverage  probability  of  a confidence  interval  on  in  terms  of  tite  parameters  N,  p, 
and  the  variant#  components  through  p = From  this  empirical  relationship 

we  can  acquire  a deeper  insight  into  the  effects  of  these  parameters  on  the  coverage 
probability  of  a confidence  interval. 

Tb  demonstrate  tiie  development  of  such  an  empirical  model,  let  us  consider,  as  an 
example,  that  h = 7.  Thus,  w (m.nj,  • • . ,117}.  We  consider  a 3^  factorial  arrange< 

meat  for  W,  d,  and  p with  the  selected  levels,  N = 50, 100,500,  d = 0.30,0.05,0.95, 
and  p = 0.1, 0.0, 0,9.  For  each  assignment  of  the  triple  (W,d,p),  designs  arc  genet' 
need  on  the  basis  of  Khuri's  (1996)  method.  Table  4.13  lists  the  number  of  designs 
generated  for  each  of  the  27  combination  of  the  levels  of  tV,  d.  and  p.  A total  of  300 
designs  were  thus  obtained.  Some  of  these  deigns  arc  listed  in  Tbble  4.19,  along  with 
the  actual  value  of  d for  a d-generaled  design. 

Four  confidence  interval  formulas,  as  described  in  Section  4.4.2,  arc  used  in  this 
investigation.  To  obtain  an  approximate  coverage  probabibty  of  a confidence  interval 


Tabic  4.19.  Generated  docigne  Tor  the  one-way  model  with  k=7  denotes  the  actual 
\iiiue  of  di  for  a ^generated  design] 


for  a gL\‘cn  design,  and  a p value,  a simulation  study  is  performed.  In  this  study, 
a = O.OS  is  used.  The  tabled  F-values  thus  used  In  the  simulation  are  follows: 
Fs7i(6,«j)  = 0.2062  Fm(6.  no)  = 2.4082 

Fs7s(43,oo)  = 0.6093  Fms(43,oo)  = 1.4890 
F«7s(93,oo)  = 0.7233  Fra(93,oo)  = 1.3224 
F,7s(493,eo)  = 1.0  Fo2s(493,oo)  = 1.0 

F«7s(6,  43)  = 0.1993  Fo„(6, 43)  = 2.7630 

F!7s(6,93)  = 0.2029  F„k(6, 93)  = 2.5658 

Ft75(6, 493)  = 0.2062  F^(B,  493)  = 2.4082 
These  values  are  obtained  by  extrapolating  the  F-tahlod  values  Burdick  and 
Graybill  (1992).  An  outline  of  the  simulation  is  described  as  follows: 

(a)  Without  any  loss  of  generality  we  can  assume  that  )i  = 0 in  model  (3.1).  We 
then  liave 

Pij  = <>i  + <u.  i = l. 2, •••,*;  j = l,2,-".7i,.  (4.49) 

Also,  without  loss  of  generality,  we  assume  4-  <7,  = 1.  Therefore,  by  recalling 

that  p=  jrfp.  Qt  ~ A'(0,p)  and  cy  ~ A'(0, 1-p). 

(b)  I^r  a given  p value  and  a given  design,  G = (nt.nj,  • • • ,n?),  a random  vector 
y is  generated  on  the  basis  of  modei  (4.49).  For  this  siinulscion,  a FORTRAN 
program  is  developed  udng  Park  and  Miller's  (1988)  uniform  random  number 
generator. 

(c)  W'ith  a simulated  data  vector  y,  two-sided  95%  con6dence  intervals  are  com- 

puted on  using  each  of  the  four  formulas  for  a cunlidenee  interval  on  [or 
equivalently,  p).  namely,  MLS,  TH,  BE,  and  MHM. 

(d)  Steps  (b)  and  (c)  are  repeated  a sufficient  number  of  times.  An  estimated  cover- 

age probability,  p,  of  a conlidence  interval  on  o’  is  thus  obtained  by  calculating 
the  proportion  of  the  intervals  that  contain  the  p value. 
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(«)  Repoar  steps  (b),  (c),  and  (d)  for  different  apecificatfon  of  a p value  and  a design 
D. 

A lota]  of  ID^DODsets  of  y in  step  (b]  were  simulated  for  each  iiietiiod.  Estimated 
coverage  probabilities  of  the  four  methods,  namely,  pui-s>  Pth>  pas,  Pmhm>  are  given 
in  Table  4.2D.  Since  several  designs  are  generated  for  each  selection  of  N and  d, 
values  of  Pi,  ^=^fLS,  TH,  BE,  and  MHM,  corresponding  to  the  same  triple  (iV,d,p) 
represent  replicated  response  values  with  N,  d.  and  p acting  as  control  variables.  The 
data  in  Thble  4.20  can  now  be  used  to  model  pi  against  N,  d,  and  p for  each  method. 
Giveo  the  nature  of  the  response  pi,  being  a probability,  it  would  be  appropriate 

consider  tbe  transformation 

a,~  = , m > 1,  ( = MLS,  TH,  BE,  and  MHM,  (4.50) 

which  maps  the  interval  (O.I)  onto  (Q.oo).  Since  several  replications  are  available 
on  P(,  and  hence  on  ie[",  t=MLS,  TH,  BE.  and  MHM,  for  each  triple  (N,^p)  in 
Table 4.19,  the  mean  iill^  and  sample  variance  ^ of  tuf*  for  triple,  i = 1, 2,*  • • , 27, 
can  be  computed.  lb  detect  a possible  linear  relationahip  between  the  mean  and 
variance  of  or  between  the  mean  and  standard  deviation  of  tv",  a simple  linear 
regression  model  without  an  intercept  using  iuTJ  and  or  tuj^  and  s"  is  fitted.  On 
the  basis  of  the  result  in  Table  4.20,  the  coefficients  of  determioalion,  B’,  for  Che 
regression  niotleLs  are  follows:  For  the  MLS  method,  the  IP  value  based  on  using 
and  sjf  U 0.66  for  m = 1,  0.81  for  m = 2 and  0-83  for  m = 3,  whereas  using  u?i 
and  s"  results  in  an  value  of  0.80  for  m = I,  0.87  for  m = 2 and  0.91  for  m s 3. 

and  0.88  for  m = 3,  whereas  using  and  s"  gives  an  B*  value  of  0.92  for  m = 1, 
0.94  for  m s 2 and  0.93  for  m s 3.  For  the  BE  method,  tbe  IP  value  for  and 
is  0.58  for  m s 1,  0.77  for  m s 2 and  0.84  for  m = 3,  whereas  using  and  s" 
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l^ble  4.2(1.  The  estimated  coverage  probability  of  o’  by  dcmlatlon  using  four  meth- 
ods for  each  of  the  generated  designs  in  Table  4.19  (nominal  confidence  coafliclenC  is 
0.95) 


products  an  value  of  0.86  for  m = 1.  0.93  for  m = 2 and  0.96  for  m = 3.  Finally, 
for  the  MHM  method,  the  fP  value  correspondinf  to  and  afT*  ^ m * 1, 

0.63  for  m = 2 and  0.74  for  m = 3,  whereas  for  5!^  and  the  iP  value  is  0.81  for 
m = I,  0.86  for  m = 2 and  0.89  for  m s 3,  This  indicates  that  there  exists  a strong 
linear  relationship  between  the  mean  and  standard  deviation  of  for  l=MLS,  TH, 
BE,  and  MH.M.  With  this  transfonnation,  a generalized  linear  model  can  be  fitted 
using  a gamma  distribution  for  itf,  IwMLS,  TH,  BE,  and  MHM.  Tb  relate  the  linear 
predictor,  4,  to  the  mean  of  w^,  a logarithmic  link  function  is  chosen.  This  link 

function  was  preferred  over  the  canonical  link  function  for  the  gamma  distribution 
because  it  produced  a more  feasible  estimate  of  the  coverage  probability.  Hence. 
4 = log(p,^),  1=MLS,  TH,  BE,  and  MHM.  Since  0<«<land0<p<l,  weuse 
the  variable  p = ^ instead  of  using  N itself  in  order  Co  reduce  the  range  of  N from 
450  to  1.  Hence,  the  region  of  interest  determined  ky  the  selected  values  of  S,  0,  and 
p in  Table  4.18  is 

C = {(P,*,fl)  I ^ < P < 0,30  < d < 0.95.  0.1  < p < 0.9}.  (4-51) 

The  linear  predictor  cotiKdered  is  4 = where  x = (p,0,p)‘,  ^ is  a vector  of 

unknown  parametera,  and  fix)  is  a vector  function  of  x.  Among  several  model  forms 
of  the  linear  predictor,  a model  having  all  two-faclor  interactions  and  a three-factor 
interaction  is  chosen  based  on  its  scaled  deviance  and  scaled  chi-squared  eCatistic. 
The  GENMOD  procedure  in  S.AS  (1997)  is  used  to  Bt  this  model.  The  fitted  model 
which  represents  the  linear  predictor,  4.  is  obtained  as  follows. 

For  the  MLS  mcihod, 

f'(x]0  = 8.4729  - 2. 744lP  + 0.3894d-  4.8564p-i-2.652lPd 

•r  2.4S40pp-i-4.73S8dp-  2.2868Pdp  (4.52) 
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For  tho  TH  morhod, 


fix)0  = 6.3928+1.2291i'  + 2.6660(>  + 3.2822p-1.3683v(i 


- 1.8059vp-3.5118dp  + 1.9789i/i>p 


{4.53) 


For  tli«  BE  inetliod, 


f\x)0  - 17.3806  - 4.4300i/-6.9012«-9.7703p+2.924lP« 


+ 4.8081PP  + 7.7741CP  - 2.9419p4p 


(4,54) 


For  the  MHM  method, 

f'{x)0  = 5.9S72  + 1.9288i'  + 3.1353d  + 4.0436p-  2.1028fd 


Tables  4.21.  4.22,  4.23,  and  4.24  show  the  standard  errors,  and  the  aasodaied 
F'values  edthe  corresponding  test  statistics  for  each  of  four  methods,  MLS,  TH,  BE, 
and  MHM,  respectively.  Note  that  the  same  model  form  is  used  to  Dt  the  data  in 
Table  4.20,  so  chat  we  can  compare  the  behavior  of  those  methods  in  terms  of  their 
(average  probability  on  which  is  the  topic  of  the  next  seciion. 

Let  ptlx)  denote  a predicted  value  ofpr,  E=MLS,  TH,  BE,  and  MHM,  at  a given  x 
inside  the  region  of  experimentation  Cbt  (4.51).  Using  models  (4.52)  through  (4.55), 
pt{x)  can  be  obtained,  at  a ^ven  point  z,  by 


for  each  method  of  f=MLS,  TH,  BE,  and  MHM,  respectively,  hlotc  that  there  arc 
several  replications  of  pt,  fsMLS,  TH,  BE,  and  MHM,  at  each  triple  (p,^,p).  There- 
fore, the  maximum  difference,  as  well  as  the  average  difference,  between  the  esti- 


- 2.7557fp-4.3187«p-h2.9445w4p 


(4.55) 


pM  = 


l+e-jfWld' 


(4.56) 


mated  coverage  probabilities  on 


Pi,  within  the  replicates  and  tiie  corresponding  predictetl  coverage  prolwhility  ob- 
tained by  model  (4.56),  which  is  ^(x).  tssMLS,  TH,  BE,  and  MHM,  at  x = (P,^.d) 
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are  calculated.  These  values  will  determine  how  good  the  Bcted  model  is  to  repre< 
sent  an  estimated  coverage  probability  on  tr^  by  simulation,  and  hopefully  the  true 
coverage  probability  on  a^.  The  predicted  coverage  probability  by  the  model,  A(z), 
and  Its  ituudmum  as  well  as  average  differences  From  p^  at  x s (p,  d,  p)  are  shown  in 
Tables  4.25,  4.26,  4.27,  and  4.28,  for  each  method  of  t=MLS,  TH,  BE,  and  MHM, 
respectively. 

We  note  that  values  of  the  maximuin  difference  range  from  -0.0348  to  Q.G447  for 
/=MLS;  from  -0.0123  to  0.0222  for  !=TH;  from  -0.0098  to  0-0179  for  (=BE;  and  from 
-0.0202  to  0.0377  for  twMHM,  whereas  values  of  the  average  difference  range  from 
-0.0102  to  0.0164  for  (=MLS;  from  -0.0078  to  0.0159  for  l=TH:  from  -0.0065  to  0.0109 
for  r=BE;  and  from  -0.0070  to  0.0302  fur  These  values  and  the  goodness 

of  fit  test  statistics  in  Tables  4.21  throu^  4.24  indicate  that  model  (4.56)  with 
IwMLS,  TH,  BE,  and  MHM,  given  by  (4.52)  through  (4.55)  represents  the  data  In 
Table  4.20  vary  well. 

Tb  see  the  effect  of  the  explanatory  variables  x on  the  coverage  probability  oil  (7^, 
contour  plots  of  Pt{x)  for  each  method  of  t=MLS,  TH,  BE,  and  MHM,  are  made. 
CoDteur  piots  for  IsMLS  are  shown  in  Figures  4.8,  4.9,  and  4.10.  Those  for  1=TH 
are  in  Figures  4.11, 4.12,  and  4.13.  Figures  4.14, 4.15,  and  4.16  are  for  l=BE,  whereas 
Figures  4-17,  4.18,  and  4.19  are  for  f=MHM. 

For  the  MLS  method,  the  predicted  coverege  probability  approaches  the  nominal 
value,  namely  0.95,  mainly  as  4i  increases  when  N is  Rxed  [see  Figure  4.8|.  For  large 
N,  d is  more  inOuetitial  than  p [see  Figure  4.8(c)],  Wlien  ^ is  Rxed,  a higher  level 
of  coverage  probability  is  obtained  by  decreasing  both  N and  p simultaneously  [see 
Figure  4.9).  If  a design  is  extremely  unbalanced  (d  = 0.30),  the  MLS  method  produces 
a smaller  coverage  probability  than  the  nominal  value,  regardless  of  the  else  of  N 
and  p [see  Figure  4.9(a)].  However,  lor  a acorly  balanced  design  (d  — 0,95),  a 
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Figure  4.9(t)l.  When  p ie  small,  by  increasing  * or  decreasing  N we  can  get  a good 
emerage  probability  [see  Figure  4.10(a)].  However,  even  when  a design  is  being  more 
unbalanced,  a desirable  level  of  coverage  probability  can  be  attaioed  by  decreasing 
the  size  of  N.  More  specificaily,  when  pis  as  small  as  0.1,  to  get  at  least  0.94  coverage 
prcbabilily  on  oj,  we  need  a design  having  ip  = 0.85  when  u = 1.0  (i.e.,  N = 450). 
If  the  design  is  unbalanced  with  t = 0.55,  it  can  provide  at  least  0-04  coverage 
probability  on  sj  when  v = 0.2  (i.e.,  N = 90).  When  p is  large,  chan^ng  N is  not 
os  effective  as  increasing  i [see  Figure  4.10(c)].  It  is  interesting  to  note  that  even 
when  fl  is  as  large  as  0.9,  the  MLS  method  gives  fairly  high  coverage  probability  on 
oj,  which  is  0.94,  for  a nearly  balanced  design  (0  = 0.9)  regardless  of  the  else  of  N. 
This  shows  that  the  use  of  the  MLS  formula  provides  a good  coverage  probability  on 
pQ  act  only  with  email  values  of  p,  but  also  with  large  valiiee  of  p os  long  as  a design 
is  nearly  balanced.  However,  across  the  range  of  the  parameter  space  we  considered, 
the  MLS  method  yields  liberal  confidence  iolervals,  particularly  when  p is  small  and 
p is  large. 

For  the  TH  method,  the  coverage  probability  approaches  tlie  nominal  value  as  p 
and  p values  are  increased  together  when  W is  fixed  [see  Figure  4.11).  This  holds 
across  all  N values  coasidered.  For  a given  Rvalue,  the  coverage  probability  increases 
as  p increases,  however,  increasing  JV  has  a little  effect  [see  Figure  4.12].  For  a 
nearly  balanced  design  {i  = 0.95),  the  TH  formuia  produces  a fairly  good  coverage 
probability  regardless  of  the  size  of  N and  p [see  Figure  4.12(c)].  In  Figure  4.12(c), 
only  one  contour  line  is  drawn  since  the  fitted  model  (4.56)  with  the  TH  formula 
produces  coverage  probabilities  close  to  0.951  when  d = 0.95  for  all  the  values  of  N 
and  0 considered.  When  the  p values  are  small  to  moderate,  the  coverage  probability 
increases  as  d increases  [see  Figure  4.15(a.b)].  FVom  Figure  4.13(a),  it  is  intsresting 
to  note  that  even  when  p is  small,  we  can  have  a good  coverage  probability  if  a dc^ 
is  nearly  balanced  (d  « 0.93)  regardless  of  the  size  of  N,  which  implies  that  the 


TH  method  U useful  for  ^ < 0.25  os  long  as  the  dralgn  is  nearly  balant'ed.  On  the 
contrary,  aoth  a large  p value,  the  coverage  probability  tends  Co  be  overestimated  for 
small  values  of  iV  and  d [see  Figure  4.13(c)].  An  approximate  coobdence  interval  with 
the  TH  method  tends  Co  be  liberal  when  both  d and  p values  are  small,  while  it  tends 
to  be  slightly  conservative  for  a large  p value. 

For  the  BE  method,  we  can  attain  the  nominal  value  of  coverage  probability  by 
increasing  p when  either  IV  or  ^ is  lixed  [see  Figures  4.14  and  4.15];  and  by  increasing 
N and*  together  when  pis  fixed  [see  Figure  4.16).  However,  the  BE  method  generally 
pnidures  from  moderatly  to  highly  conservative  confidence  intervals  except  wlien  p 
is  large  [see  Figure  4.16(c)). 

The  behavior  of  the  MHM  confidence  interval  is  very  similar  to  that  of  TH’s  [sec 
Figures  4.17,  4.18,  and  4.19). 

In  summaty,  the  MLS  method  is  good  for  large  values  of  * with  sinnll  p values, 
however,  it  produces  n highly  liberal  coofidence  Inlorval  when  * is  smalJ,  or  p is  large 
unless  * is  large.  Both  TH  nnd  MHM  methods  generally  produce  a good  coverage 
probabiltty  unless  both  * and  p are  smali  and  N is  small  to  moderate.  In  this  case, 
however,  MLS  is  better  than  TH  and  MHM.  The  BE  method  is  good  only  when  p is 
large,  or  both  N and  * are  large  if  p is  small. 

4.4.4  Comparison  of  methods  in  terms  of  coverate  orobalalitv 

The  MLS  procedure  works  well  if  is  small.  However,  if  is  large,  this  procedure 
may  result  in  very  liberal  intervals  (Burdick  and  GraybiU,  1992,  p.  70).  Although 
Burdick  and  Eickman  (1986)  indicated  chat  this  procedure  also  yields  liberal  intervals 
in  extremely  imbalanced  demgns,  the  combined  effect  of  the  variance  components  and 
the  demgn  on  this  procedure  has  not  yet  been  studied.  Tbomas  and  Hultquist  (1978. 
p.  5S6)  stated  that  “they  would  not  expect  exact  methods  to  give  mudi  better  results, 
at  least  for  the  designs  considered’' , that  is,  better  results  tban  their  formula  gave. 


Figure  4.9.  Contour  plots  of  the  predicted  coverage  probability  on  with  model 
14.52)  for  each  level  of  v = ^ using  the  Modified  Urge  Sample  (MLS)  method 
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Figure  4i(.  Contour  plots  of  the  predicted  coverage  probability  on  oj  with  model 
(4.53)  for  each  level  of  d (t'=^l  using  the  Modlbcd  Large  Sample  (MLS)  method 


Figure  4.11.  Contour  plotrs  of  the  predicted  coverage  probability  on  a\  with  model 
(4.53)  for  each  level  of  r-  = using  the  Tlioniae-Hultqiijst’a  (TH)  method 


Figure  4.12.  Contour  plots  of  the  predicted  coverage  probability  on  with  model 
(4.53)  Tor  each  level  of  * (i^^)  using  the  Thomaa-Hultquist’a  (TH)  method 
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FiguK  4.13.  Contour  plots  of  the  predicted  coverage  probability  on  with  model 
(4.53)  for  each  level  of  p [i^®]  using  the  Thomas-HuJlquist's  (TH)  method 


Figure  -i-H-  Contour  plots  of  the  predicted  coverage  probability  on  a\  with  model 
(4.&4)  for  each  level  of  e = using  the  Burdick- Eickman'a  (BE)  metbod 
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Figure  4.15.  Contour  plots  of  tbe  predicted  coverage  probability  on  rr^  with  model 
(4.54)  fur  each  level  uf  4 [t>=^|  using  the  Biirdick-Eicknuui’e  (BE)  method 


Figure  4.16.  Contour  plots  of  the  predicted  uoverage  probability  on  with  model 
(4.54)  for  each  level  of  p |v=^|  using  the  Burdick-Eidunao's  (BB)  method 
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Figure  4.17.  Contour  ploto  of  tbo  predicted  coverage  probability  on  al  with  madel 
(4.05)  for  each  level  of  i/  s ^ uabg  the  Modified  Harmonic  Mean  (MHM)  method 
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Figure  4.1S.  Contour  plots  of  the  predicted  coverage  probability  on  with  model 
(4.65)  for  each  level  of  d using  the  Modihed  Harmonic  Mean  (MHM)  method 


Figure  4-19.  Coutour  plots  of  tbe  predicted  cut'etage  probability  on  with  model 
(4-bo)  for  each  level  of  p using  the  Modified  Harmonic  Mean  (MHM)  method 
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since  the  appronmalion  of  the  H-'-raiidoin  varisble  with  the  chi-sqiiared  distribution 
worked  very  well  for  most  one-way  unbalanced  designs  they  had  comndered.  Burdiclc 
end  Elclunan  (1986)  showed,  using  simulation  study,  that  their  procedure  provided 
a higher  coverage  probability  on  than  ThoioaS'Kultquist's  procedure  across  the 
experimental  situations  they  considered.  They  indicated  that  their  procedure  could 
always  be  recommended  over  the  Thomas-Hultquisl’s  method  (Burdick  and  Eickman, 
1986,  p.  218).  They  showed  tliat  their  interval  had  a coverage  probability  that  was 
generally  at  least  as  large  as  the  nominal  value.  It  was  also  noted  that,  although 
interval  (4.47)  was  more  conservative  than  Thomas-Hultquist's,  these  two  methods 
produced  confidence  intervals  whose  lengths  did  not  di^r  by  mote  than  5%.  Khuri 
(1999)  showed  that,  by  replacing  oa  with  run,  the  approximation  associated  with  the 
IV-random  variable  in  Thomas  and  Hultquist  (1978)  is  adequate  when  oj  is  large. 
It  was  also  shown  that,  when  oj  is  small,  the  approximation  with  n„,  was  slightly 
better  than  that  with  for  upper  toil  probabilities,  while  the  former  was  slightly 
worse  than  the  latter  for  lower  tail  probabilities.  Ke  also  compared  the  use  of  ns 
versus  that  of  ru.  with  regard  to  the  power  of  the  F-leat  concerning  oj. 

interval  on  o^  using  their  estbnoted  coverage  probabUilies  obtained  by  eimulation. 
Partial  comparisons  based  on  contour  plots  using  predicted  coverage  probabilities 
obtained  by  the  6tted  models  were  already  mode  in  the  previous  section.  Table  4.29 
shows  the  rsnges  of  simulated  coverage  probabilities  resulting  from  the  replicates 
for  a given  triple  (fV,0,p).  Note  that  the  nominal  confidence  coefiiclent  is  0.95. 
Using  the  nonnal  approximation  Co  the  bbiomial,  if  Che  true  confidence  coefficient  is 
0.95,  there  is  less  than  2.5%  chance  that  an  estimated  confidence  coefficient  based 
on  10,000  replications  will  be  less  than  0.9457  (=0.95-1,96^^^^^^).  Also,  there 
is  Ires  than  24%  chance  that  an  estimated  confidence  coefficient  will  be  more  than 
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0.9543  (=0-95+1.96^^*  Value*  which  are  significantly  below  the  slated  level 

in  each  triple  {N,  d.  p)  are  indicated  using  a symbol  *.  and  tJioee  with  ngnificantly 
above  the  stated  level  ate  iudicated  using  n symbol  #.  The  MLS  procedure  generaUy 
piovides  liberal  confidence  intervals,  but  is  good  when  d is  large  and  p is  small.  It 

iV  is  not  large.  The  TH  and  MHM  procedures  behave  similarly,  and  produce  liberal 
intervals  when  p is  small.  When  p is  moderate  to  large,  the  TH  procedure  teaches 
the  nominal  value  of  coverage  probability  unless  N ie  small-  However,  the  MHM 
procedure  tends  to  produce  a slightly  higher  coverage  probsbility  than  the  nominai 
value  for  the  moderate  to  large  values  of  p,  especially  when  N is  large.  The  BE 
ptoiodurc  provides  too  conservative  intervals  when  p is  small.  This  procedure  is  good 
with  moderate  to  large  values  of  p if  d Is  not  too  small. 

Pearson's  chi-square  goodness-of-fit  statistic  is  evaluated  over  300  generated  de- 
signs, namely, 

y (p,.-0.9sy  ^ y 

in  order  to  see  the  agreement  of  the  estimated  coverage  probability,  obtained  by 
simulation,  with  the  nominal  value  over  the  re^on  of  experimenution  C in  (4.51). 
The  values  of  are  = 0.5156,  X^  = 0.0224,  Xgg  = 0.06561,  and  = 
0.0356.  The  Xj  statistics  evaluated  Just  over  the  replicates  for  each  triple  (W,d.p) 
are  listed  in  Thble  4.30. 

The  average  value  of  the  estimated  coverage  probability  obtained  by  simulation 
within  the  replicates  for  each  triple  {N,  i.  p)  is  supposed  to  represent  the  true  coverage 
prubabilitv  for  the  same  triple.  Plots  of  the  average  estimated  coverage  probability 
against  N for  each  combination  of  d and  p values  are  given  in  Figure  4.20.  In  these 

plots,  the  dashed  line  ( ) is  for  the  MLS  method,  the  dotted  Une  { ) is  for 

the  TH  method,  the  short  broken  line  ( — ] is  for  tbe  MHM  method,  and  the  loog 


in&l  conR- 
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broken  line  ( ) is  for  the  BE  method.  The  MLS  method  piodiices  closer  coverage 

prohsbUity  to  the  nominal  level  than  other  methods  with  small  to  moderate  fiise  of 
•V  when  p is  small.  However,  if  S is  large,  the  TH  and  MHM  methods  are  better 
than  Ilic  MLS  method  when  p Is  small.  Fhr  moderate  to  large  values  of  P,  the  TH, 
BE.  and  MHM  procedures  provide  good  coverage  probabilities  regardless  of  the  size 
of  .V  and  ®,  Figure  4.21  is  for  the  plots  of  the  average  coverage  probability  against  ® 
for  each  combination  of  jV  and  p values.  A better  performance  of  the  TH  procedure 
than  the  MHM  procedure  can  be  seen  when  JV  = 50  and  p = O-l.  or  when  Af  = 100 
and  p s C.l,  for  all  values  of  ®.  In  both  cases,  however,  the  MLS  method  produces 
better  coverage  probability  than  both  TH  and  MHM.  When  fV  = 500  and  p = 0,1,  the 
MHM  method  is  best,  though  it  is  slightly  better  than  the  TH  method.  The  coverage 
probability  of  the  MLS  procedure  increases  toward  the  nominal  value  as  ® increases 
regardless  of  the  siee  of  N and  p,  however,  iu  performance  is  poor  compared  u>  the 
other  procedures.  For  large  values  of  p,  no  distinction  can  be  made  among  the  TH, 
BE  and  MHM  procedures.  Finally,  plots  of  the  average  coverage  probability  against 
p for  given  values  of  N and  ® are  given  in  Figure  4.22.  The  coverage  probability  of 
the  MLS  and  BE  decrease  as  p increases.  The  TH  and  MHM  methods  behave  in 
a mrailar  pattern.  They  perforin  lietter  than  the  BE  method  in  general.  They  also 
produce  better  coverage  probability  than  the  MLS  method  except  for  soialJ  p and 
small  to  moderate  values  of  JV  regardless  of  ®. 

Based  on  the  above  plots  we  can  summarize  the  behavior  of  the  four  methods  as 
follows:  The  MLS  method  is  preferred  over  the  other  methods  when  p is  small  and 
N is  small  to  moderate.  Otherwise,  the  TH  and  MHM  methods  give  better  coverage 
probability,  and  the  distini'tion  iietween  these  two  methods  is  small.  The  BE  method 
is  also  ss  good  as  tlie  TH  and  MHM,  hut  only  when  p is  large. 
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Figure  4.20.  The  average  estimated  coverage  probability  within  the  replicates  against 
V (v  = for  each  combinatian  of  levels  of  4I  and  p.  Lince  repreaant  the  MLS 
method  ( ),  the  TH  method  ( — ),  the  MHM  method  { — ),  and  the  BE  method 
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Fig^ire  4.21.  The  average  eatimated  coverage  probability  witbin  the  repbeatee  againat 
« for  each  combination  of  leveia  of  v and  p (v  = ^).  Lines  represent  the  MLS  method 
i ),  the  TH  method  (•••*),  the  MHM  method  ( ),  and  the  BE  method  ( ~) 


Concludine  remarks 


The  coverage  probabiJity  of  a confidence  interval  depends  not  only  on  the  design 
used,  but  also  on  the  values  of  the  variance  components.  Several  apprccdmate  con- 
fidence intervals  on  for  the  one-way  random  mtxlei  are  available  in  the  statistical 
literature.  However,  their  performance  in  terms  of  the  coverage  probability  as  it  re- 
lates u>  the  variance  components  and  the  design  has  not  been  studied  in  the  statistical 

In  tills  study,  vm  Investigate  the  elTect  of  imbalance  as  well  as  the  effect  of  variance 
components  on  the  coverage  probability  of  four  approximate  confidence  intervals.  To 
do  this,  tor  each  of  the  confidence  intervals,  an  estimated  coverage  probability  is  ob- 
tained by  simulation  with  a given  design  and  a specified  ratio  of  variance  components. 
Thus,  this  probability  is  a function  of  the  total  number  of  observations,  the  degree  of 
imbalance,  and  the  ratio  of  tbe  variance  components  for  a fixed  number  of  treatment 
levels.  This  probability  is  empirically  morleled  in  terms  of  these  vaiiablK  using  gener- 
alized linear  models  techniques  for  eacli  method.  The  fitted  model  is  used  to  observe 
the  efiect  of  these  variables  on  the  coverage  probability.  The  four  methods  can  then 
be  compared  on  the  basis  of  their  respective  models.  Although  no  »nglc  method  is 
considered  best,  the  modified  large  sample  interval  is  useful  when  p is  small  with 
nearlv  balanced  dt«ign.  It  also  performs  better  than  tbe  other  intervals  when  p is 
small  and  N is  small  to  moderate.  The  Thomas-Hultqulst  and  tbe  modified  harmonic 
mean  intervals  perform  well  for  moderate  to  large  values  of  p regardles,  of  the  maes  of 
jV  and  i.  The  Burdick-Eiekman  interval  is  also  good  for  large  values  of  p,  however, 
it  is  too  conservative  if  p is  small. 


CHAPTER  5 
FURTHER  RESEARCH 


Three  extensions  of  the  use  of  the  QDGs  in  Khurt  (1997)  were  given  in  Chapter  3, 
namely,  the  applications  to  the  one-way,  two-way  crttss-cl&seification  without  interac- 
tion, and  thrce-foid  nested  classification  random  models;  the  use  of  unbalanced  data, 
and  the  ron^eration  of  maximum  likelihood  (ML)  estimation.  Fbrther  exten^ons 
ate  needed.  For  example,  consideration  of  other  methods  of  estimating  variance  com- 
ponento  such  as  REML  and  truncated  ANOVA.  The  application  of  the  methodology 
with  mixed  effects  models  are  also  of  interest.  Moreover,  consideration  of  models 
with  random  effects  having  distributions  other  than  the  normai  is  of  interest  and  can 
provide  added  flexibility  to  the  method  of  comparing  designs  on  the  basis  of  QDGs 
or  EQDGs. 

fn  Chapter  -1,  we  investigated  the  effect  of  Imbalance  on  various  quantities  of  in- 
terest. They  were  the  probability  of  getting  a negative  ANOVA  esuraate  of  trj,  the 
power  of  the  F-teat  for  er^  with  liomogcncity  assumption  concerning  ths  error  van- 
ances,  the  power  of  the  same  test  under  heterogeneity  assumption  of  error  variances, 
and  the  coverage  probability  of  a confidence  interval  on  These  Investigations  were 
carried  out  using  an  unbalanced  one-way  random  model.  An  extension  of  the  pro- 
posed methodology  to  higher-order  models  is  of  interest.  Extending  Khuri’s  (1937) 
method  of  generating  designs  for  the  one-way  model  to  crossed  elaasification  models 
is  straightforward.  For  the  nested  cla^cation  model,  however,  the  extension  can 
be  more  involved  since  different  degrees  of  imbalance  exist  at  different  stages  in  the 
nesting  relationship.  Although  designs  can  be  generated  separately  for  each  stage, 
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we  need  to  have  an  overall  meaeure  of  imbalance  for  all  the  stages  in  order  to  see  the 
effect  of  imbalance  of  a design  on  a quantity  of  interest.  The  effect  of  imbalance  on 
other  qnantitia  is  also  of  interest.  Fbr  example,  the  small  sample  behavior  of  ML 
estimator  can  be  ereamined  by  modeling  the  mean  squared  error  of  the  estimator  in 
terms  of  the  variance  components  and  the  design-  With  a nested  claasificntion  model, 
an  investlgatioD  of  the  so-calied  "counter-balancing  effect”  in  Cummings  and  Gaylor 
(1P74),  which  was  brieliy  reviewed  in  Section  2.1,  is  ofinlereet. 


APPENDIX 

SPLVS  CODE  FOR  GENERATING  DESIGNS 
Following  is  the  SPUJS  program  lo  generate  a design  for  the  given  value  of  the 

treatment  group  la  6xed  by  t = 5 in  this  program.  For  other  values  of  k,  the  func- 
tion, "imbal.gct.d'',  should  be  modified. 


att(»U<M*StS(tHtSiaU  program  begins  *M«tn«*SttlH*H«*aSMSS«« 


Q[i.J3  <-  lMrt(J«(j*D) 
0Ci.lt]  <-  l/sqrtCk) 

> 

Qfk.ld  <-  l/aqrttk) 
forCl  in  S;kH 
3 <-  l-l 

DCl.j]  <-  -j/aqrt<ivj) 
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'ormal  randoo  numbers 
onlfonaly 
'or  the  algoritha 


[ to  generate  r: 
generating  thle  randoo  ' 


out  <-  natrljt(-999,lt,l) 
r <-  8qrtCU“phi)/(hephi)> 

flag  <-  0 

ubileCflageeO  t iteration<lOO(X)0){ 
iteration  <-  iteration  • 1 
H <-  raonB(kl,0,l) 

■ <-  Bqrt(t(«)X»WI) 
for<i  in  i:kl){ 

Z[iJ  <-  r*S[i]/e 

) 

Z[k]  <-  l/eqrt(k) 

Q <-  imbal .heloertCk) 

J <-  Q X<X  Z 
if<oin«)>0){ 
flag  <-  flagtl 


Calculate  all  possible  coobinatione  of  2‘k  values  on  vector  d tpei) 
Thle  function  is  specifically  for  k«5  C2‘k^2  vaya) 

Dbal.get.d  <-  fuoctionCk2,k,bcH 
JCk  5) 

atap('*you  chose  an  incorrect  function.  This  is  Juet  for  k*6“) 

forCJl  in  i;2H 
for(j2  in  1;2H 
for(j3  in  1:2)< 
fer(J4  la  l;2H 


chktflag.J  <-  tordi 
flagl  <-  Jlag-1 
lfCflagl—l>){ 
outCl,:  <-  cMod 


ford  in  I:fl»gl){ 

dupld.llagJ  <-  6«a(al)B(clii[flag,]-cUiCl,])) 

> 

if  tany{lopl[.fla*3  ■■(»)•( 

flag  <-  flag-1 


program  anda 


tyminuatn  aampla  run  uul  tha  follouiag  paraiators:  #»#»tMaw«» 

«»»««»»*»»•  It'S,  l!*50.  phl-0.8.  rliO"0.6,  and  nanf?  »*«»*»«#>•# 


•atmtMtMttfttWMMMtM 
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L2.J 

C3.1 

»,] 

[S.l 

[6,] 

[7,] 

[8,: 

C9.3 

[10,] 


C.I]  [,2]  [.3]  [,4]  [,B] 


C.a  [,7]  C,6]  [.8] 

14  2 8 0.793SS0B 

1 13  14  0.7987220 

17  8 13  0.7987220 

11  11  3 0.8012821 

6 12  7 0.7987220 

15  16  4 0.7987220 

3 17  9 0.8012821 

a 10  9 0.8038685 

10  4 10  0.7987220 

12  3 IS  0.7961783 


tr  the  actnal  phi-value  for  a phi'senaratail  desiga 


Remarks: 

(1)  in  order  to  generate  a design  for  an  unbalanced  one-way  model,  It  is  only  needed 

to  epetdfy  jV  end  d values.  A p value  does  not  used  to  generate  a dosign.  Tbc 
p value  is  used  when  we  ralrulate  the  values  oT  the  quantities  of  Interest  wbirdi 
are  considered  in  this  dissertation. 

(2)  When  this  program  seems  not  ninning  properly  or  it  takes  too  long  to  wait,  I 

recommend  to  iocreese  the  sire  of  “d”,  or  to  decrease  the  else  of  “want",  or  to 
increase  the  size  of  “S'. 
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